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A  Fast  Algorithm  of  the  Finite  Difference  Method  for  the 
Computation  of  the  Transonic  Flow  Past  an  Arbitrary  Airfoil 
with  the  Conservation  of  the  Full  Potential  Equation 

Huang  Mingke 

(Nanjing  Aeronautical  Institute) 

I.  Introduction 

There  are  two  primary  finite  difference  methods  to  solve  the 
full  potential  equation  for  the  transonic  flow  around  an  airfoil. 

One  of  them  was  developed  by  Garabedian,  Korn  and  Jameson^ ^-2]^ 

This  method  first  conformally  maps  the  exterior  of  the  airfoil  on 
the  physical  plane  onto  the  interior  of  a  circle.  Then,  a  fast 
direct  method  and  a  linear  relaxation  technique  are  used  in  the 
circle  to  solve  the  full  potential  finite  difference  equation. 
Because  the  center  of  the  circle  corresponds  to  infinity  on  the 
physical  plane,  the  velocity  potential  is  infinite  at  that  point. 

In  order  to  avoid  this  singularity,  Jameson  introduced  the 
perturbation  potential  to  eliminate  the  free  flow  velocity 
potential  which  includes  such  singularity  to  complicate  the 

[3] 

governing  equation.  The  other  method  was  developed  by  Holst 
He  used  an  arbitrary  curve  coordinate  which  matches  the  surface 
to  transform  the  computation  zone  into  a  rectangle.  Then,  the  AF2 
iteration  scheme  is  used  to  solve  the  finite  difference  equation 
in  the  rectangular  region.  AF2  is  a  fast  and  effective  iteration 
scheme  which  has  already  been  widely  used^^”^^.  Because  Holst 
used  an  arbitrary  non-orthogonal  curved  coordinate,  it  is  highly 


1 


flexible.  However,  it  takes  some  computer  time  to  create  such  a 
mesh.  In  addition,  the  equation  becomes  more  complex  after  the 
transformation. 

On  the  basis  of  these  two  methods,  a  simpler  method  is 
developed  in  this  work.  We  first  employ  conformal  mapping  to 
transform  the  computing  region  into  a  unit  circle.  Then,  a 
radial  transformation  is  carried  out  in  the  circle  to  reduce  the 
effect  of  singularity  of  velocity  potential  at  the  center  of  the 
circle.  Hence,  it  is  not  necessary  to  introduce  the  perturbation 
potential.  The  governing  function  is  thus  much  simpler  than 
those  used  by  Jameson  and  Holst.  Then,  the  equation  of 
conservation  of  full  potential  is  discretized  using  a  center 
difference  scheme.  Moreover,  an  artificial  density  is  introduced 
to  the  artificial  viscosity  in  the  supersonic  region^"*  ’  . 

Finally,  the  finite  difference  equation  is  solved  by  the  AF2 
scheme  developed  by  Holst.  A  computer  program  has  already  been 
written  for  this  algorithm.  Computations  have  already  proven 
that  this  algorithm  is  successful  and  effective. 


II.  Computation  Method 

First,  the  exterior  profile  of  the  airfoil  on  the  complex 

physical  z-plane  is  mapped  onto  the  unit  circle  a  by  conformal 

[21  [81 

mapping.  The  derivative  of  the  transformation  function  is  ’ 


dz 

da 


•/* 


(1) 
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This  paper  was  received  on  August  1,  1983  and  revised  manuscript 
was  received  on  October  19. 

“  i  0 

where  e  is  the  trailing  edge  angle  of  the  airfoil.  Let  a  =  e  / 20 

/ r  and  0  =  0  corresponds  to  the  trailing  edge  of  the  airfoil. 

This  transformation  is  able  to  map  an  open  trailing  edge  onto  an 

enclosed  circle.  The  coefficient  cn  in  equation  (1)  needs  to  be 

[81 

determined  by  iteration  ,  usually  by  7-8  iterations.  Because 
the  fast  Fourier  transform  (FFT)  technique  is  used  in  the 
iteration,  the  conformal  mapping  process  takes  less  than  a  half 
minute  on  a  TQ-16  computer.  Then,  it  is  transformed  according  to 
£=1/o  to  map  from  outside  onto  the  inside  of  the  unit  circle. 
Therefore,  r,  0  and  £  are  the  polar  coordinates  on  the  plane. 

The  velocity  and  velocity  components  are  rendered 
undimensional  by  using  the  critical  sonic  speed.  The  density  is 
made  dimensionless  with  the  value  at  the  stationary  point.  The 
conservative  full  potential  equation  on  the  £  plane  thus  becomes 


where  <t>  is  the  velocity  potential.  The  velocity  components  U 
and  V  along  the  0  and  r  coordinates  on  the  physical  plane  can  be 
expressed  as 

,,  r  d<t>  v  gl_-  dt 

L  =  H  dr  /  is 


where  H  =  |dz/dc|  is  the  transformation  mode.  The  relation 
between  p  and  the  velocity  components  is: 


(L''rP)|’" 
y  + 1  J 


(4) 
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where  ir  is  the  specific  heat.  The  boundary  conditions  on  the  S 
plane  are 


Along  the  surface  ( 3^/ar)r_-j=0 
At  infinity  or  when  r-  0, 

_,C0Svg  j 6 — tan'1  [(l-A/i)1'’  tan c6~a  —  ^)2 

where  V  and  M  are  the  incident  velocity  and  Mach  number  of  the 
flow,  o  is  the  angle  of  attack  and  u  is  the  imaginary  part  of 
the  coefficient  c„  ,  which  corresponds  to  the  zero  attack  angle  of 
the  airfoil  at  low  velocities.  This  can  ensure  that  the 
conformal  mapping  is  defined.  r  is  the  ring  around  the  airfoil, 
which  is  determined  by  the  Kutta  condition.  The  Kutta  condition 
requires  that  (a^/ae)Q_Q  =  0- 

From  equation  (6)  we  find  that  d~(1/r)  near  the  origin. 
Therefore,  the  boundary  conditions  cannot  be  directly  written  for 
r=0.  However,  it  is  possible  to  choose  the  external  boundary  on 

[Q1 

the  physical  plane  at  a  far  away  location  and  the  flow 
parameters  there  are  considered  as  the  free  flow  values.  If  the 
external  boundary  conditions  are  rewritten  on  the  5  plane  at  r=r„ 
where  rmin  is  very  small,  because  the  radial  step  Ar  in  the 
finite  difference  method  may  be  of  the  same  order  of  magnitude  as 
rmin»  t^ie  finite  difference  of  approximation  of  the  derivative 
ad/dr  has  a  very  large  error  near  rmin.  For  this  reason,  Jameson 
introduced  a  perturbation  potential  to  overcome  this  difficulty. 
We  employ  the  following  radial  transformation 


(5) 

(6) 


The  value  at  each  radial  mesh  point  can  be  calculated  in  the 
advanced.  Therefore,  it  does  not  add  any  complexity  to  the 
governing  equation. 

Figure  2  shows  the  physical  plane  and  various  transformation 
planes.  In  the  physical  plane,  a  seam  from  the  trailing  edge  of 
the  airfoil  is  dividing  the  computation  region  into  singularly 
connecting  regions  to  made  the  velocity  potential  singular  in  it. 
The  jump  of  4>  across  the  seam  is  equal  to  the  circular  moment. 

The  computation  region  in  the  (R,e)  plane  is  rectangular.  In 
order  to  create  the  difference  lattice,  it  is  divided  into  imax 
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meshes  in  the  ©-direction  and  j  meshes  in  the  R-direction.  F 

max 

corresponds  to  Hence,  the  computation  region  in  the  (R, 

9)plane  corresponds  to  0<i<i  +1  and  j  .  „ ,  as  shown  in 

max  Jmin  J  Jmax 

Figure  2(c). 


Figure  1.  The  Radial  Coordinate  Transformation  Function 
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Figure  2.  Physical  Plane,  Transf ormation  Plane  and  Computation 
Plane 


1 .  external  boundary 

2.  surface 

3.  surface 

4.  external  boundary 

5.  trailing  edge 

6.  surface 

7.  upper  seam 

8.  lower  seam 

9.  external  boundary 

(a)  physical  plane  z=x+iy 

(b)  transformation  plane  £ 

(c)  computation  plane  (R,e) 

[3] 

In  analogy  to  the  Holst  method  ,  the  finite  difference 
approximation  of  equation  (8)  is 

U>.,  =  3  ,  [(-rf-nr)  , 

}(R)>  ■  J  (10) 

where  6Q,  6^  and  6Q,  6^  are  the  forward  and  backward  difference 
operators.  Because  the  curved  coordinate  for  conformal  mapping 
is  a  coordinate  to  fit  the  surface  of  the  object,  the  streamline 
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direction  in  the  supersonic  region  is  relatively  close  to  the  0- 
direction.  Therefore,  we  only  introduced  the  artificial 


viscosity  in  the  0-direction.  To  this  end,  an  artificial 
density  p  is  used  in  equation  (10)  to  replace  p.  We  get 

p  .  -  v,.j.  i  ( 1 


where  v  =  max[0 , (M^-1 )c]  (12) 

Here,  M  is  the  local  Mach  number.  The  coefficient  c  is  between 
1-2.  Another  limit  is  that  the  value  of  v  cannot  be  larger  than 
1.  In  equation  (11)  k  is  chosen  to  be 

k  =  j  0  uhen  ut^,j>0 
1  when  Ul+%(j<0 

When  we  calculate  the  density  p  on  a  mesh  point  based  on 

equation  (5),  the  required  U  and  V  are  calculated  by  the  center 

difference.  The  densities  with  the  1/2  indicator  were  computed 

by  taking  the  average.  Due  to  the  fact  that  p=p  in  the  subsonic 

region,  therefore,  the  difference  approximation  equation  (10)  has 

a  second  order  accuracy.  However,  the  difference  term  with 

respect  to  0  only  has  a  first  order  accuracy. 

r  3  ] 

The  iteration  scheme  can  be  written  as 

NCn+wUn=0  (13) 

where  n  is  the  number  of  iterations,  Cn=tfn+"'-tfn  is  a  correction 
term,  L^n  is  the  residual  value,  and  u>  is  the  relaxation  factor. 
The  operation  N  should  be  selected  to  be  as  close  to  L  as 
possible.  It  should  also  facilitate  the  process  to  find  its 


m 


inverse.  We  chose  a  two  multiplying  factor  format 


aA'C;.,=  —  [a—  J,(p/(i?)) 


where  a  is  the  acceleration  parameter  (not  to  be  confused  with 
the  front  attack  angle). 

Each  iteration  is  specifically  implemented  in  two  steps: 
First  step: 


Second  step: 


c<z-  d 

(14) 

1 

(15) 

where  f?  .  is  an  intermediate  result.  In  the  first  step,  it  is 
1  >  J 

only  necessary  to  solve  two  diagonal  linear  algebraic  equation 

sets  for  each  e=const.  line.  To  this  end,  it  is  also  required  to 

specify  the  boundary  conditions  for  f?  . .  Similar  to  that  in 

1  *  J 

reference  [3],  it  is  assumed  that  af/3R=0  along  the  surface. 

With  the  intermediate  value  f?  . ,  it  is  possible  to  solve  the 

1  >  J 

tridiagonal  algebraic  equations  for  each  R=const.  line  from  the 

external  boundary  inward  in  the  second  step.  Rigorously,  this 

equation  is  not  a  tridiagonal.  There  is  one  non-zero  element  on 

the  upper  right  side  and  one  on  the  lower  left  side^^. 

Therefore,  the  velocity  potential  in  the  n+l^  iteration  is  ^n+^  = 

— • 

<t>n  +  Cn  to  complete  an  iteration.  In  the  second  step,  +oB60is  the 

term  to  provide  convergence  to  the  term  in  the  supersonic 

ro  .nn[2] 
region 


The  initial  field  for  iteration  is  chosen  as 


^  =  Vm  (V-r-p-  jcos(<2~0-*i)-i-^-  tan'1  [(1  -  Mi)1"  tan  (0  -  a  -  m>3 

where  the  first  term  on  the  right  side  is  the  velocity  potential 

of  an  incompressible  circulationless  flow.  The  second  term  is 

the  circulation  part.  In  this  equation,  a  again  represents  the 

angle  of  attack.  Obviously,  with  the  exception  that  the 

circulation  is  not  yet  determined,  this  initial  field  has  already 

satisfied  the  boundary  conditions  at  the  external  boundary  as 

well  as  those  on  the  surface  of  the  object.  The  initial 
0 

circulation  r  is  chosen  to  make  the  initial  field  meet  the  Kutta 
condition.  When  the  attack  angle  is  not  too  large,  we  get 

r « -  r  =  i!tF-  sip 
~  (1  (l-Mi)*'* 

For  an  incompressible  flow,  the  initial  field  chosen  is  the  final 
solution.  For  a  subsonic  flow,  the  initial  circulation  is  close 
to  the  final  value. 

Because  the  initial  field  already  satisfies  the  Kutta 
condition  and  the  circulationless  part  of  the  initial  field  also 
meets  the  external  boundary  condition,  therefore,  the  correction 
field  obtained  in  each  iteration  is  zero  at  the  external 
boundary.  In  order  to  determine  the  correction  for  circulation, 
the  correction  field  is  required  to  satisfy  the  Kutta  condition. 
Let  us  assume  that  we  have  already  found  the  following 
circulations  rn,rn”^,  rn~^ .  in  earlier  iterations.  We 


n+  <1 

estimated  rn  in  the  north  iteration  by  extrapolation.  The 
increment  of  circulation  is  estimated  to  be 

Ar,=<arf2r,-3r,*‘+r*:,> 


where  «r  is  the  low  relaxation  factor.  This  circulation 
increment  can  be  used  as  an  auxiliary  relation  in  solving 
equation  (15): 


C;.,=c-  :  at-,,  c-  ci.,— &r. 


After  solving  C?  .  from  equation  (15), the  circulation  of  the 
1 » J 

correction  field  is 


A  r  =  c* 


Therefore,  rn+^=rn+Ar,  which  is  used  to  correct  the  circulation 
on  the  external  boundary  for  the  next  iteration. 


III.  Numerical  Examples 

On  the  (R,0)  plane,  the  computation  region  is  divided  into 

20  meshes  along  the  R -direction,  i.e.,  j _ =40  and  j  .  =20.  In 

the  6-direction,  imax  Is  equal  to  64  and  128.  The  acceleration 
parameter  chosen  is  a  geometric  sequence  formed  by  8  numeric 
values.  In  each  iteration,  one  of  the  values  in  the  sequence  is 
used  in  order.  After  a  complete  cycle,  it  starts  all  over  again 
cij  is  used  to  represent  the  minimum  in  the  sequence  and  <*n  is 
the  maximum.  When  a^  is  small,  it  converges  rapidly.  If  it  is 
too  small,  the  correction  in  each  iteration  is  too  large  to  lead 
to  the  solution  to  diverge.  Therefore,  a^  has  a  lower  limit. 
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The  selection  of  the  optimal  acceleration  parameter  must  rely  on 
trial  computation.  In  the  following  example,  let  us  choose  <*^  =  10  ~ 
20  and  the  ratio  of  two  neighboring  numbers  in  the  geometric 
series  a  is  1.6.  In  the  initial  stage,  a  smaller  a^  is  chosen  to 
accelerate  the  development  of  the  flow  field.  In  the  later 
stage,  a^  is  chosen  to  be  larger  to  avoid  the  solution 
diverging.  Other  parameters  are  chosen  as  follows:  &=0'N-1.0,  w 
=  1.8,  <ur=0 . 7~1 . 0 ,  c=1 . 

Figure  3  shows  the  calculated  pressure  distribution  on  the 

NACA  0012  airfoil  when  M^O.63  and  attack  angle  a=2°.  This 

example  belongs  to  the  subcritical  case.  The  mesh  chosen  is  64 

x  20.  The  result  is  in  agreement  with  that  obtained  by  using  the 

T  3 1 

Holst  method  with  a  149  x  28  mesh  .  This  example  needs  123 
iterations  to  obtain  the  convergent  result  shown  in • the  figure. 

In  this  case,  the  pressure  coefficient  does  not  vary  by  more  than 
0.001  over  8  iterations,  requiring  approximately  35  minutes  on 
the  computer . 

Figure  4  shows  the  results  obtained  with  both  coarse  and 
fine  meshes  for  the  same  airfoil  when  M  =0.75  and  a=2*.  This  is 

w 

an  example  of  the  supercritical  case,  which  requires  277 

iterations.  When  a  128  x  20  mesh  is  chosen,  the  result  is  in 

[33 

total  agreement  with  that  of  Holst's 
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Figure  3.  Calculated  Results  of  NACA  0012  Airfoil  at  Mm=0.63 
and  a  =  2°  Using  a  Coarse  Mesh  (64x20) 
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Figure  4.  Calculated  Pressure  Distribution  on  NACA  0012  Airfoil 
at  M  =0.75  and  a=2° 


(a)  coarse  mesh  (64x20) 

(b)  fine  mesh  (128x20) 

1 .  this  work 

2.  chord  coordinate 

3.  this  work 

4.  chord  coordinate 

All  the  above  computations  were  done  on  Chinese  made  709  and  / 24 
TQ-16  computers.  It  should  be  pointed  out  that  these 
computations  were  made  by  choosing  the  acceleration  parameter  and 
other  parameters  conservatively  because  of  insufficient  computer 
time.  Preliminary  computation  already  demonstrated  that  these 
parameters  have  an  extremely  large  effect  on  the  rate  of 
convergence.  Therefore,  there  is  a  large  potential  to  accelerate 
the  process. 
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FLOW  PAST  AN  ARBITRARY  AIRFOIL  WITH  THE 
CONSERVATIVE  FULL-POTENTIAL  EQUATION 

Huang  Mingke 
Wanting  Aeronautical  Institute > 

Abstract 

Based  on  the  methods  developed  by  Jameson  and  Holst,  a  computer 
program  has  been  developed  for  computation  of  the  transonic  flow  past  an 
arbitrary  airfoil  by  the  finite  difference  method.  A  conformal  mapping  is 
applied  to  map  the  exterior  of  the  airfoil  onto  the  interior  of  a  circle. 
By  a  radial  transformation  reducing  the  effects  of  the  singularity  at  the 
centre  of  the  circle,  the  use  of  the  perturbation  velocity  potential  is 
avoided.  The  governing  equation  simpler  than  those  used  by  Jameson  and 
Holst  is  approximated  by  a  finite  difference  equation,  which  is  then  solved 
by  AF2  iteration  scheme  in  computing  plane.  The  computations  of  the 
pressure  distribution  over  the  airfoil  of  NACA  0012  for  subcritical  and 
supercritical  cases  show  the  results  in  excellent  agreement  with  those,  by 
Holst's  method. 
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A  Mixed  Finite  Difference  Analysis  of  the  Internal  and 
External  Transonic  Flow  Fields  of  Center  Cone  Inlet 
Luo  Shijun,  Shen  Huili,  Ji  Minggang  and  Xing  Zongwen 
(Northwestern  Polytechnical  University) 

and 

Dong  Songye  and  Han  Aiqing 
(31th  Institute,  Ministry  of  Astronautical  Industry) 

I.  Introduction 

In  reference  [1],  the  internal  and  external  transonic  flow 

field  of  a  Pitot  tube  type  axisymmetric  inlet  was  successfully 

calculated  by  applying  a  mixed  finite  difference  method  to  the 

transonic  steady  potential  of  the  flow.  The  shape  of  a  center 

cone  gas  inlet  mostly  belongs  to  the  contraction-expansion  type. 

The  disturbance  on  the  internal  wall  surface  is  relatively 

intense.  In  addition,  due  to  the  presence  of  the  centerbody,  the 

transversal  relaxation  line  which  is  able  to  transmit  outlet 

resistance  against  pressure  disturbance  is  cut  off.  Therefore, 

the  alternating  direction  linear  relaxation  iteration  technique 

used  in  the  calculation  of  a  Pitot  tube  inlet  could  not  be 

[  1  ] 

successfully  applied  .  In  this  work,  on  the  basis  of 
experience,  measures  such  as  multiple  regional  iteration, 
choosing  low  relaxation  factor  in  the  initial  stage,  and  giving 
the  appropriate  initial  field,  were  taken  to  obtain  stable  and 
convergent  results.  Specifically,  the  method  used  in  this  work 
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was  based  on  a  large  longitudinal  disturbance  (small  transverse 
disturbance)  flow  potential  equation.  A  mixed  difference  scheme 
was  used.  A  directional  derivative  method  was  used  to  accurately 
insert  the  surface  boundary  conditions.  We  either  started  with  a 
zero  initial  field  or  with  a  known  solution  to  a  similar 
condition,  or  with  a  solution  of  a  one-dimensional  entropy  flow 
as  the  initial  internal  flow  field  to  solve  the  difference 
equation  by  using  a  regional  multi-layer  relaxation  iteration 
method.  Finally,  the  precise  potential  flow  was  used  to  make 
corrections  (i.e.,  to  use  the  small  transverse  disturbance 
potential  solution  as  the  initial  field  to  solve  the  precise 
velocity  potential  equation) .  This  method  was  used  to  calculate 
the  inlet  described  in  reference  [2]  as  well  as  inlets  of  other 
geometric  shapes.  The  incident  flow  Mach  number  is  M<j>=0 . 72~1 . 27. 
The  outlet  Mach  number  is  M2=0 . 205~0 . 82 .  The  calculated  median 
flow  velocity,  pressure  coefficient  distribution,  and  the  shape 
and  position  of  shock  wave  are  in  agreement  with  experimental 
results . 


II.  Basic  Equations 

When  the  Mach  number  of  the  incident  flow  is  not  too  large 
(e.g.,  less  than  1.3),  the  iso-entropic  inviscid  assumption  is 
valid  for  a  steady  flow  in  an  axisymmetric  inlet.  The  velocity 
potential  equation  is: 
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where  x  and  y  are  the  axial  and  radial  direction  of  the 

cylindrical  coordinate,  respectively,  u  and  v  are  the  axial  and 

radial  components  of  the  flow,  respectively,  a  is  the  local 

sonic  speed.  <t>  ,  <t>  ,  d  ,  dvv  and  <i>  are  the  first  and  second 
x  y  xy  xx  yy 

order  derivatives  of  the  disturbed  velocity  potential  function  0 
(x,y)  with  respect  to  x,  y,  respectively.  The  local  sonic  speed 
a  is : 

This  paper  was  received  on  August  26,  1983. 

+  +  (2)  / 

In  equation  (1),  p  is  a  coefficient.  When  p=1,  equation  (1) 
is  the  precise  velocity  potential  equation.  When  p=0 ,  equation 
(1)  becomes  a  large  disturbance  velocity  potential  equation  in 
the  x-axis  direction. 

ao  represents  the  sonic  speed  of  the  incident  flow  in  the 
expression  for  a.  qo  is  the  incident  flow  velocity  and  y  is  the 
adiabatic  index. 

III.  Division  of  Computation  Region  and  Difference  Lattice 
Because  the  difference  computation  is  calculated  on  a 
physical  plane,  the  size  of  the  region  to  compute  will  affect  the 
convergence  and  stability  of  the  results.  When  the  incident  flow 
is  subsonic,  shock  waves  are  only  formed  locally  near  the  cowl. 
The  difference  computation  is  relatively  stable  and  the  computed 
area  needs  not  be  too  laree.  When  the  incident  flow  is 


supersonic,  the  front  far  field  boundary  only  needs  to  include 


the  front  of  the  separated  bow  wave.  However,  the  side  field  and 
rear  field  must  be  sufficiently  large  to  ensure  the  stability  of 
the  separated  bow  wave.  Furthermore,  it  should  not  hinder  the 
acceleration  of  the  supersonic  flow  behind  the  wave  in  the  inlet 
to  become  supersonic.  In  this  case,  the  far  field  boundary  might 
be  selected  to  be  around  1000R,  where  R  is  the  cross-section  of 
the  inlet  (see  Figure  1).  When  the  computed  area  is  too  small, 
it  may  lead  to  the  oscillation  or  even  divergence  of  the 
solution . 

The  lattice  is  divided  orthogonally  along  the  x  and  y 
direction.  In  order  to  make  the  lattice  points  also  full  on  the 
body  surface,  the  spacings  are  not  equal  (see  Figure  2).  Near 
the  leading  edge  of  the  inlet,  the  meshes  are  closer.  The  mesh 
spacing  is  usually  determined  according  to  a  law  of  algebraic 
series.  The  ratio  of  two  neighboring  lattice  spacing  is 
generally  controlled  at  under  2.  Local  points  may  exceed  this 
limit. 


♦  *  0 


Figure  1.  Computed  Region  and  Boundary  Conditions 


Figure  2.  Unequal  Lattice  Spacing 


In  order  to  reduce  the  computing  time,  a  method  to  gradually 
decrease  the  spacing  of  the  meshes  is  used.  Initially,  a 
relatively  loose  mesh  (42x51)  is  used.  After  the  flow  field 
stabilized,  the  meshes  are  tightened  (81x100)  in  order  to  obtain 
better  flow  field  parameters. 

IV.  Boundary  Conditions 

1.  Far  Field  Boundary  Condition.  The  disturbed  velocity 
potential  in  the  ir  front  field,  side  far  field,  and  rear  far 
field  outside  the  inlet  is  given  to  be  ^=0. 

The  given  flow  velocity  in  the  field  far  away  from  the 
outlet  is  q2* 
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2.  Centerline  Boundary  Condition.  Because  of  axisymmetry, 
v=«Sy=0  at  the  centerline. 

3.  Centercone  Leading-edge  Boundary  Condition.  The 
stationary  point  condition  is  used  for  all  points  on  the  leading- 
edge  of  the  centercone,  i.e.,  u=^>x  +  qa)=0,  v=0 . 

4.  There  are  two  ways  to  treat  the  boundary  conditions  for 
the  leading-edge  of  the  cowl. 

(1)  u=^x+qoo=0,  based  on  a  blunt  leading-edge. 

(2)  Based  on  a  sharp  leading-edge,  the  point  is  split 
into  two.  The  upper  point  is  located  on  the  upper  surface  and  is 
treated  as  a  body  boundary  point  of  the  upper  surface.  It  is 
considered  as  an  extension  of  the  upper  surface  in  the 
computation.  Furthermore,  the  <t>  value  on  the  extended  surface  is 
ext:  ipolated.  The  lower  point  is  located  on  the  lower  surface 
and  is  treated  as  a  body  surface  boundary  point  of  the  lower 
surface.  In  the  computation,  it  is  assumed  to  be  an  extension  of 
the  lower  surface.  Furthermore,  the  <b  value  on  the  extended 
surface  is  extrapolated. 

5.  Body  Surface  Boundary  Conditions.  Let  us  assume  that 
the  equation  of  the  body  surface  is  y=f(x).  A  directional 
derivative  method  is  used  to  accurately  meet  the  body  surface 
boundary  condition,  i.e., 

rfy-f'Cx)  <q„+*x>  (3) 

If  s  is  used  to  express  the  direction  along  the  body  surface,  the 
directional  derivative  equation  can  be  written  as 

<tfg  =^xcos(x,s)  +  <£gcos(y,s)  (4) 
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where 


COS  ( X ,  s)  =  -  - 

v/  1-f'U)1 
cos  (y,  s)  = — -JL.S*} - 

v'  1  -f'(x)1 

By  summarizing  the  above,  we  get 


d,  =  _ 


<t>, 


v''  ”  1-r  /'(*)' 

<t>,=  (q (x) 


(5) 

(6) 


The  second  order  derivative  of  the  perturbed  velocity  potential 
on  the  body  surface  can  also  be  expressed  in  terms  of  directional 
derivatives  : 


4,.  =  <*..v'Y- /'(*}’  -*„/'(*> 

®„=  j£,,  -i,./'  (*>]%/  1-r/'  (*>‘  (■*)’ 


(7) 

(8) 


where  the  second  derivative  is  expressed  as  the  following 

using  an  analytical  extension  method: 


Oye 


2  [  _±\  (for  upper  surface)  (9) 

'  Ay,  V  Ay,  V 

^  y,'‘  (for  lower  surface)  (10) 


Here,  the  subscripts  i  and  j  represent  the  lattice  numbers  of 
body  surface  points. 

In  these  equations,  the  formulas  to  calculate  the  partial 

derivatives  <b  ,  <t>  ,  and  4>  of  perturbation  velocity  potential 

s  x  s  y  s 


In  the  direction  tangent  to  the  body  surface  are  as  follows 
(using  the  positive  slope  surface  as  an  example) 


d,  -  As?  -  As? 

As. As,., (As, -As,.,) 

AS,  As,., 


(11) 


(12) 


'  As,- As,., 


(13) 


^  if ( As,.,  +  As,-,) 1  As* -  —  ^,-i./-i(As,-|-4~Asj-,) l+^i.,,i.|Ai! ., 

Asi-,As,_a(As,_i  +  As,*—]) 

JL  —  «• ,  ~  (^i) 

Vfl  , 

As,., 


(14) 


(15) 


where 


i  —  ^»)  i>  i  “*  C^f)  i-l. 
V»i  “  "x  . 

As,., 

Aj.^n/  Ax’,  4- A j/j 


(16) 


(17) 


_As,_,«  >/  A*J_,+  AvJ-i 
As,.,  -  %/  Ax’-,+  Aj/< 


(18) 


(19) 


Equations  (11),  (12)  and  (13)  are  for  subsonic  points,  and  (14), 
(15)  and  (16)  are  for  supersonic  points. 


V.  The  Mixed  Difference  Scheme 
Different  difference  schemes  for  x-direction  partial 
derivatives  are  used  for  supersonic  and  subsonic  points  in  the 
flow  field.  For  a  supersonic  point,  the  difference  of  the 


x-direction  partial  derivative  employs  a  downwind  scheme.  For  a 
subsonic  point,  a  center  scheme  is  used.  The  y-direction  partial 
derivatives  employ  the  center  difference  scheme.  With  the 
exception  of  some  special  cases ,  a  second  order  accuracy  scheme 
is  used  for  all  first  derivatives  and  a  first  order  accuracy 
scheme  is  used  for  all  second  derivatives  in  order  to  match  the 
scheme  to  the  accuracy  desired.  The  main  formulas  are: 


.  _  ^..,[(Ax^,+  A*,.,)*- Ax;_,]-^,_l>,(Ax,_l  + 

1  Ax.^Ax,.,  (Ax.^  +  Ax,.,) 

i  _«  <bi.,  A*,.,-<A,.i.i(A*,-i  +  Ax,.,)  A*,., 

"  Ax,_,  Axj_,(Ax,.,  +  Ax,_,) 


(20) 

(21) 


4  _ ± - 

(Ay,+  Ay*.,)  (AxUl+  Ax,.t)  Axi.Axi., 

x((4i,i+,-4i.j-i)C<Axi.l  +  AxI_,)*-Ajr,i.J 

(0i-l,  (.1 I -l)  (Ax,_,+  Ax,.,)1 

+  (^i-J./  +  l  Axj  } 

J  ^  Ax?.,  +  ^,.,-  (Ax?-Ax?_,) -^,.,.,Ax» 

*  Ax,.,  Ax,  (Axi_,  +  Ax.) 


(22) 


(23) 


^  » 2  Ax,.,-^..,  (Ax,.,  +  Ax,)+0,.,.,Ax, 

Ax,_,Axp  (Ax,_,  +  Ax,) 


(24) 


(Ay<  +  Ay/.,)  (Ax,_,+  Axf)  Ax,_,  Ax( 

^  C  i+i  — ^  AXj  - 1  ^  (fit  i+i  /— 

f  -  Ayf.,+^,.,  (Ay?~Ay?.,)-^,.)-tAy? 

4y*  Ay*.,  (Ay,+  Ay,.,) 


(25) 


(26) 


a  (Ayj_,+ Ayj)+^„t-i  Ay< 

”  Ay<+,  Ay(  (Ay,.,+ Ay<) 


(27) 
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where  equations  (20),  (21)  and  (22)  are  in  the  downwind  mode  and  / 2 9 
(23)-(27)  are  in  the  center  mode. 

VI.  Velocity  Determination  and  Treatment  of  Shock  Wave  and  Sonic 
Speed  Points 

In  establishing  the  difference  equation  for  a  point  in  the  flow  field, 

it  is  necessary  to  determine  the  velocity  at  that  point  to 
see  whether  it  is  subsonic  or  supersonic  in  order  to  employ 
a  comparable  difference  scheme  in  iteration.  The  determination 
process  is  shown  as  follows: 


shock  wave  supersonic  supersonic  sonic 

point  point  point  point 


where  the  subscript  s  indicates  that  a  one-sided  formula  is  used 
and  the  subscript  0  indicates  that  the  center  formula  is  used. 


Sonic  and  shock  wave  points  are  treated  as  subsonic  points. 


For  a  sonic  point,  notice  that  equation  (1)  becomes 


d  +  d  / y=0 

y  y  ^  y  '  J 


(28) 


VII.  Multi-regional  Multi-layer  Linear  Relaxation  Iteration 


Method 

If  the  control  function  (1)  is  written  as 

R.d  +  Rod  +  Rod  +  R/d  =0  (29) 

1  xx  2*yy  3  xy  4  y 

We  get  the  following  by  substituting  equations  (20)— (27)  into 


(29). 


Ad •  .  !  +  Bd,  ,  +  Cd.  , . !  =  D  (30) 

If  equation  (30)  is  written  for  every  point  on  all  vertical 
lattice  lines  (i=constant)  and  the  coefficients  A,  B,  C  and  D  are 
known,  then  a  series  of  tridiagonal  equations  is  the  result. 
Hence,  it  is  possible  to  solve  this  problem  by  a  chasing  method. 

In  a  center-body  conical  inlet,  the  internal  is  mostly  a 
contracting-expanding  conduit.  With  a  supersonic  incident  flow, 
when  the  computation  begins  with  a  zero  initial  field,  it  is 
equivalent  to  calculating  the  flow  field  from  a  supersonic  flow 
field.  In  this  case,  the  wall  disturbance  in  the  expanding  tube 
will  lead  to  the  further  acceleration  of  the  flow.  Thus, 
calculated  results  will  further  deviate  from  the  actual  flow. 
However,  the  boundary  condition  at  the  outlet  requires  the 
presence  of  subsonic  flow  in  the  expanding  channel.  For  this 
reason,  when  the  entire  field  is  linearly  relaxed  along  the  y- 
direction,  the  wall  disturbance  inside  the  tube  and  the 
disturbance  of  the  reverse  pressure  at  the  exit  are  seriously  out 


of  coordination.  In  response  to  this  special  characteristic, 
wall  disturbance  must  be  suppressed  in  the  initial  stage  of  the 
calculation  in  order  to  accelerate  the  propagation  of  the  inverse 
pressure.  However,  because  the  center-body  blocked  the 
disturbance  of  inverse  pressure  propagating  along  the  x-direction 
relaxation  line,  when  we  used  the  alternate  direction  linear 
relaxation  iteration  method  to  calculate  the  experimental  data, 
it  did  not  converge  over  long  iterations.  For  this  reason,  a 
regional  multi-layer  lower  order  relaxation  iteration  method  was 
used  in  this  work,  i.e.,  the  entire  flow  field  was  linearly 
relaxed  in  the  y-direction.  In  addition,  lower  relaxation  was 
used  in  the  initial  stage.  The  entire  flow  field  was  divided 
into  internal  and  external  flow  areas.  Or,  the  internal  flow 
could  be  further  divided  into  a  contracting  section,  an  expanding 
section  and  a  straight-section.  Moreover,  the  number  of 
relaxation  iterations  for  the  internal  flow  might  be  increased 
(e.g.,  one  iteration  for  external  flow,  1-5  iterations  for 
internal  flow  and  1-5  times  more  iterations  for  the  expanding 
section).  Until  the  flow  field  is  essentially  stable,  then  the 
entire  flow  field  was  treated  as  a  whole  by  relaxation  iteration. 

Our  experience  indicates  that  this  method  not  only  can  ensure  the  /30 
stability  of  the  computation  but  also  can  accelerate  convergence. 

This  is  because  it  not  only  facilitates  the  propagation  of  the 
inverse  pressure  disturbance  at  the  outlet  but  also  can  adjust 
the  coordination  between  wall  disturbance  and  outlet  pressure 
disturbance.  Consequently,  an  initial  flow  field,  similar  to  the 
actual  flow  field  can  be  established.  This  facilitates  the 
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solution  of  a  flow  field  starting  from  a  zero  initial  field 
because  its  structure  is  not  known. 


Conventional  or  improved  iterations  may  be  used  in  the 
linear  relaxation  iteration  method.  The  improved  iteration  not 
only  can  speed  up  the  rate  of  convergence  but  also  will  not  cause 
a  supersonic  flow  field  to  become  divergent.  Different 
relaxation  factors  were  used  in  the  supersonic  and  subsonic 
regions.  A  relaxation  factor  w  is  defined  as  follows: 


where  the  superscript  (n)  represents  the  number  of  iterations  and 

4)*.^^  represents  the  result  of  n*"*1  iteration  not  yet  treated  by 
1  *  J 

relaxation . 


VII.  Initial  Field  and  Convergence  Standard 
The  selection  of  the  initial  field  is  very  important.  It 
will  directly  affect  the  convergence  and  time  of  computation.  It 
is  feasible  to  use  a  zero  initial  field  (i.e.,  the  disturbance 
velocity  potential  ^=0)  for  the  external  flow  of  the  inlet  due  to 
the  fact  that  it  is  relatively  simple.  It  also  makes  the 
computation  relatively  more  stable.  However,  for  the  internal 
flow  field  of  the  inlet,  if  a  one-dimensional  flow  method  can  be 
used  to  estimate  an  initial  flow  field,  then  the  rate  of 
convergence  can  be  increased.  If  this  is  difficult  to 
accomplish,  then  we  can  also  start  from  a  zero  initial  field.  In 
this  case,  because  the  initial  field  is  far  different  from  the 


final  solution,  relatively  smaller  relaxation  factors  must  be 
used  over  a  longer  period  of  time  to  perform  iterations  in  layers 
during  the  initial  stage  of  computation  to  prevent  it  from 
becoming  divergent. 

Whether  an  iteration  is  convergent  must  be  evaluated  by  the 
computing  process.  Let  us  assume  that  the  following  is  true  for 
all  points 

r,:\  U 

Then,  the  standard  convergence  can  be  specified  as: 

A*  saO-’-lCT’ 

In  this  work,  in  addition  to  the  above  convergence  standard, 
the  main  standard  of  convergence  used  is  that  the  relative  error 

of  any  cross-sectional  flow  rate  is  less  than  e  (e=0.0'K0.05) : 

A 

where  A4^n^  =  I  <|^ n^-  I  ,  is  the  gas  flow  rate  across  the 

1  i  e  'max  i 

i^  cross-section  and  4>e  is  the  flow  rate  across  the  exit  cross- 
section.  Our  experience  shows  that  it  is  a  convergence  criterion 
which  reflects  the  physical  nature  of  the  flow.  It  not  only  can 
avoid  a  false  convergence  due  to  using  excessively  small 
relaxation  factors  but  also  can  prevent  any  misjudgement  caused 
by  velocity  potential  fluctuation  due  to  the  pressure  of  sonic 
line  or  small  shock  wave  vibration  when  the  incident  flow  is  near 
sonic  or  supersonic. 


IX.  Examples 

The  method  introduced  above  was  used  to  calculate  the 
following  flow  conditions  of  the  inlet  with  a  conical  centerbody. 

(1)  Mw=1.2278,  M2=0.3481  (E) 

(2)  Mm=1 . 2278 ,  M2=0.2489  /3 

(3)  M^l.2278,  M2=0.2053 

(4)  M^l.0453,  M2=0.3435  (A) 

(5)  M^l.0453,  M2=0.3025 

(6)  M^l.0453,  M2=0.2056 

In  addition,  we  also  calculated  internal  and  external  transonic 
flow  field  of  the  first  type  (60-40)  inlet  reported  in  reference 
[23. 

The  computation  was  carried  out  on  a  SIEMENS  7760  computer 
using  two  sets  of  lattices,  i.e.,  42x51  and  81x100.  The 
criterion  of  convergence  used  is  ( A<|» / 4»e ) <  1  —  3% .  The  numbers  of 
iterations  required  to  obtain  a  convergent  solution  are 
approximately  320  and  360,  respectively,  taking  about  500  and  800 
seconds  of  CPU  time.  For  the  first  type  of  inlet  reported  in 
reference  [23,  we  get  the  following  solutions: 

Small  transverse  disturbance  (u=0)  solutions  with  42x51  and 
81x100  meshes,  as  well  as  the  accurate  velocity  potential 
solution  (y=1)  with  42x51  meshes  when  the  incident  flow  Mach 
number  M  =1.27  and  the  flow  coefficient  CA=0.655. 

09 

Small  transverse  disturbance  solution  with  42x51  meshes  when 
M  =1.27  and  CA=0.721. 

09 

Small  transverse  disturbance  solution  with  42x51  meshes  when 
M  =1.27  and  CA=0.762. 
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[2] 

Comparisons  of  major  results  to  experimental  data  are 
shown  in  Figures  3^8.  They  show  that: 


Figure  3.  Pressure  Distribution  Along  Centerbody  and  Cowl 
Surface  (M  =1.27,  M  =0.5413,  CA=0.655) 

1 .  centerbody 

2.  theoretical  cone  surface  pressure  ratio 

3.  cowl 

4.  calculated  value  (p=0,  42x51) 

5.  experimental  value  reported  in  reference  [2] 


Figure  4.  Pressure  Distribution  Along  Centerbody  and  Cowl 
Surface  (M  =1.27,  M  =0.6329  CA=0.721) 

«  on 


1 .  centerbody 

2.  theoretical  cone  surface  pressure  ratio 


4.  calculated  value  (m=0,  42x51) 

5.  experimental  data  reported  in  reference  [2] 


Figure  5.  Pressure  Distribution  Along  Centerbody  and  Cowl 

Surface  (M  =1.27,  M  =0.7058,  CA=0.762) 

00  00 

1 .  centerbody 

2.  theoretical  cone  surface  pressure  ratio 

3.  cowl 

4.  calculated  value  (m=0,  42x51) 

5.  experimental  value  reported  in  reference  [2] 


Figure  6.  Comparison  of  Pressure  Distribution  Along  Centerbody 
and  Cowl  Surface  CM  =1.27,  M  =0.5413,  CA=0.655) 
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1 .  centerbody 

2.  theoretical  cone  surface  pressure  ratio 


3. 

4. 

5. 
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meshes 

meshes 


experimental  data  reported  in  reference  [2] 


1.  For  a  supersonic  inlet,  when  the  center  conical  angle 
and  cowl  leading  edge  angle  are  not  too  large,  the  results 
obtained  by  using  a  mixed  difference  method  (such  as  pressure 
distribution  on  surface  and  position  and  shape  of  the  bow-shaped 
shock  wave)  with  a  small  transverse  disturbance  steady  transonic 
velocity  potential  agree  well  with  the  experimental  data  (See 
Figures  3~5  where  Dc  is  the  diameter  of  the  cross-section  of  the 
entrance  of  the  inlet) .  The  accuracy  is  already  acceptable  for 
engineering  applications  when  using  the  42x51  mesh  size.  There 
is  no  need  to  make  corrections  by  tightening  the  meshes  or  using 
the  full  velocity  potential  equation  (See  Figure  6). 

2.  Figure  7  shows  that  the  position  of  the  shock  wave  is 

determined  from  the  pressure  distribution  curves  (Figures  3-5) 

based  on  the  pressure  ratio.  The  calculated  value  is  1.405. 

When  the  experimental  data  was  processed,  it  was  chosen  to  be 

[21 

1.383  by  taking  loss  into  account  .  The  position  of  shock  wave 
obtained  by  this  method  is  located  one  mesh  ahead  of  the 
experimental  result.  This  is  due  to  the  use  of  a  non¬ 
conservative  difference  scheme  (see  references  [3]  and  [4]). 

3.  Near  the  leading-edge  of  the  cowl,  the  calculated 
pressure  is  higher  than  the  experimental  value.  It  is  so  even 
with  the  solution  to  the  accurate  velocity  potential.  This  is  so 
even  with  the  solution  partially  separated  near  the  leading  edge 
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under  these  three  conditions.  The  total  pressure  loss  is 
significant.  Assumptions  such  as  inviscidy  and  isoentropy  are  no 
longer  valid.  However,  the  separated  area  is  very  small  and  will 
not  affect  the  pressure  drag  of  the  entire  inlet.  The 
engineering  value  of  this  method  will  not  be  hurt. 

4.  The  calculated  intensity  of  the  conical  shock  wave  and 
the  cone  Mach  number  (pressure)  are  in  good  agreement  with 
theoretical^^  and  experimental  results  (see  Figures  3,4,5  and 
8).  Because  the  meshes  are  loose,  the  position  of  the  shock  wave 
cannot  be  easily  determined.  However,  the  error  is  less  than  one 
mesh  space. 
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Figure  7.  Distance  Between  Shock  Wave  and  Cone  Tip  vs.  Flow  Rate 

(M  =1.27) 
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1.  calculated  value 

2.  experimental  value 
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Figure  8.  Pressure  Distribution  Along  Centerbody  and  Cowl  vs. 
Flow  Rate 

1.  theoretical  cone  surface  pressure  ratio 

2.  centerbody 

3.  cowl 


X.  Conclusions 

1.  In  this  work,  we  successfully  calculated  the  internal 
and  external  flow  fields  of  a  supersonic  inlet  based  on  the  large 
disturbance  axial  velocity  potential  equation  by  using  a  simple 
non-rotational  mixed  difference  scheme,  accurate  surface  boundary 
conditions  and  the  regional  linear  relaxation  iteration 
technique.  This  method  has  good  convergence  and  stability.  The 
accuracy  and  required  computer  time  can  meet  practical 
engineering  needs. 

2.  When  calculating  the  flow  field  of  an  inlet  with 
complicated  geometry,  conventional  linear  relaxation  iteration 
techniques  cannot  easily  render  convergent  solutions  because  the 
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flow  parameters  in  various  regions  are  not  in  coordination.  This 
problem  can  be  overcome  by  using  multi-layer  iteration  and  using 
a  low  relaxation  factor  in  the  initial  stage  or  by  choosing  an 
initial  field  closer  to  the  real  solution.  The  experience  gained 
in  this  work  can  prepare  us  to  extend  our  work  to  calculate  the 
internal  and  external  flow  fields  of  complicated  three- 
dimensional  inlets. 
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A  MIXED  FINITE  DIFFERENCE  ANALYSIS  OF  THE 
INTERNAL  AND  EXTERNAL  TRANSONIC  FLOW 
FIELDS  OF  INLETS  WITH  CENTERBODY 
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Abstract  ,  ~ 

,  A  mixed  finite  difference  method  for  calculating  the  external  and 
internal  flow  field  around  inlet  with  centerbody  is  presented.  First,  calcu¬ 
lation  by  mixed  finite  difference  method  of  the  velocity  potential  equation, 
with  small  disturbance  in  the  transverse  direction  using  cartesian  mesh, 
irroiationa!  schemes  and  exact  body  surface  boundary  conditions  is  carried 
out  to  obtain  a  basic  field  solution  including  the  shape  and  location  of 
the  shock  and  the  sonic  line.  Then,  the  full  potential  equation  is  used  to 
improve  the  accuracy  of  the  computed  value  of  field  variables.  The  use  of 
multi-layer  line  relaxations  along  the  radial  lines  is  effective  for  inlet 
with  centerbody,  and  in  this  case,  more  relaxation  sweeps  are  carried  out 
(with  smaller  relaxation  factor)  inside  the  inlet  than  the  relaxation  swe¬ 
eps  carried  out  outside  the  inlet.  Computations  have  been  made  for  axisvm- 
metric  inlet  with  different  f reest ream  Mach  nu mbers  1 ,04~  1 .27.  Com¬ 
putation  results  show  that  the  method  is  promising. 
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Numerical  Method  for  Viscous  Shock-layer  Near  The  Stagnation 

Line  of  a  Blunt  Body 
Yu  Qingwen  and  Wang  Ruquan 
(The  Computer  Center,  Academia  Sinica) 

I.  Introduction 

When  studying  a  hypersonic  viscous  flow  around  a  blunt 
object,  the  prediction  of  aerodynamic  and  thermodynamic 
parameters  in  the  stagnation  region  is  also  very  important.  When 
the  viscous  layer  is  very  thin,  the  classical  boundary  layer 
equation  is  always  used.  When  the  mutual  interference  between  a 
viscous  gas  and  an  inviscid  flow  is  considered,  the  boundary 
layer  model  is  no  longer  applicable.  The  viscous  shock-layer 
equation  can  be  used  in  the  flow  field  of  a  wide  range  of 

[i  2] 

incident  Mach  number  and  Reynolds  number  ’  . 

In  recent  years,  many  people  have  dedicated  their  effort  on 
the  computation  of  the  flow  field  near  the  stagnation  line .  [  3>4,5 ] Seme 
of  them  crossed  the  shock  wave  and  some  separated  the  shock  wave. 
However,  when  the  shock  wave  was  separated,  they  neglected  the 
effect  of  the  downstream  flow  field  on  the  stagnation  region, 
which  usually  resulted  in  large  error.  Because  the  flow  near  the 
stagnation  line  is  in  the  subsonic  range  in  the  shock  layer,  the 
physical  equations  describing  the  flow  characteristics  are 
essentially  elliptical.  Hence,  when  the  effect  of  the  downstream 
flow  on  the  stagnation  line  is  considered,  it  very  often  leads  to 
new  unknown  parameters  in  the  equations.  Thus,  there  are  more 


unknowns  than  equations  to  make  the  problem  non-unique.  Just  for 
this  reason,  many  people  discarded  parameters  relevant  to  the 
downstream  flow  in  their  studies. 

In  this  work,  based  on  the  viscous  shock  layer  equation,  a 
set  of  second  order  normal  differential  equations  near  the 
stagnation  line  is  derived.  In  order  to  reflect  the  effect  of 
downstream  flow  on  the  stagnation  point,  the  shock  wave  boundary 
condition  is  expanded  into  a  Taylor  series  around  the  stagnation 
point,  including  an  unknown  parameter e 2 =(1/2) . (d  l /dx^) .  It  is 
related  to  the  curvature  of  the  segmental  shock  wave  near  the 
stagnation  point  which  remains  unknown  unless  the  flow  field  at 
the  head  is  solved.  We  used  the  curvature  of  an  ideal  shock  wave 
to  approximate  it  and  the  numerical  results  agree  with  the 
overall  solution. 


II.  Basic  Equations 

The  basic  equations  describing  the  flow  characteristics  near 
the  stagnation  point  are  the  viscous  shock  layer  equations^ ^  . 

Near  the  stagnation  point,  the  solution  can  be  expanded  into  the 

[6] 

following  • 

u  =  u.j  (y)x  (2.1) 

v  =  -( 1-%k2x2)vQ(y)  (2.2) 

p  =  ( 1-k^x^)pQ(y)-%x^P2(y)  (2.3) 

T  =  (1-k2x2)  T0(y)  (2.4) 

Manuscript  received  on  March  23,  1983,  revised  on  June  22.  This 
paper  was  presented  in  the  1978  Meeting  of  Computational  Aero¬ 
dynamics  . 
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p  =  P0(y)  (2.5) 

m  =  M()(y)  (2.6) 

A  =  A0(y)  (2.7) 

2 

e  =  eQ+e2x  + .  (2.8) 

The  last  equation  is  the  expansion  of  the  distance  of  shock 
wave  separation  at  the  stagnation  point.  A  series  of  normal 
differential  equations  cna  be  obtained  by  substituting  the  above 
equations  into  the  viscous  shock  layer  equation.  After  using 
factors  such  as  R,  Vw,  P#D,  Wo,  \m,  pjf^,  WJP/R,  and  R/W^  to 

render  length,  velocity,  density,  mean  molecular  weight,  thermal 
conductivity,  viscosity,  pressure,  temperature  and  isobaric 
specific  heat  non-dimensional,  the  viscous  shock  layer  equation 
in  the  stagnation  region  is  as  follows: 


1  d  du,  w  _  niti.  \  du,~  r  (kp.v, -/£),«,) 
Re  dy  V'  dy  )  \p‘l*  RT)  ~dy  +  [ - ~H - 


(2.9) 


.  r.  2 p-u.  A 

'TcPct,--~-  =  0 

d.tL^l.ov  -  o 

dy  '  dy 


Pi^PcT, 

pl  =  pt(7') 


(2.10) 

(2.11) 

(2.12) 

(2.13) 

(2.14) 

(2.15) 


where  H=1+ky,  fl0 =2k/ ( 1+ky) ,  6=Cp/RePr,  and  (ll,la)  are  control 
parameters.  (0,1)  is  the  thin  shock  layer  equation,  (0,0)  is  the 


/  ■  /v-v!vl  vlvl*  /vly.'  r  .y  .'r  Ay  o~.y -I 
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boundary  layer  equation  and  (1,1)  is  the  viscous  shock  layer 
equation. 

In  a  chemical  equilibrium,  it  is  necessary  to  change  the 
equations  of  energy  and  state.  In  this  case,  Cp  and  W  are  not 
constants.  Instead,  they  are  functions  of  p  and  T.  We  employed 
an  approximate  thermodynamic  function  table  method  to  calculate. 
When  the  gas  is  in  equilibrium,  the  equations  of  energy  and  state 
are  as  follows: 


a^(*^)+(*rp‘t'‘+Jr‘<W  =  °  (2.16) 


dy 


p0  =  p0  T„  /W0  (2.17) 

where  h^,  and  hp  represent  the  derivatives  of  enthalpy  with 
respect  to  temperature  and  pressure.  In  summary,  the  stagnation 
line  equation  is  a  set  of  quasi-linear  differential  equations 
with  ux  ,  v0  ,  p„  ,  pa  ,  T„  ,  po  and  p0  as  unknowns. 


III.  Boundary  Conditions 

Let  us  assume  the  following  conditions  exist  on  the  body 
surface : 


ux  =  0 
Vo  =  0 

T»  =  T, 


(3.1) 

(3.2)  / 

(3.3) 


The  shock  wave  is  considered  as  an  interrupted  surface  and 


the  classical  Rankine-Hugoniet  condition  is  the  external  boundary 
condition.  It  is  expanded  into  a  Taylor  series  at  the  stagnation 
point.  By  using  equation  (2.8),  the  expressions  on  the  shock 


wave  can  be  obtained  as  follows: 


y-1  \  M  l)l  +  ke. 


2Z-L+ _ 2 _ 

y+1  (y+l)Afi 


(3.4) 

(3.5) 


P. 


Pi. 


-/y-*  +  2  V1 

*  Vy+l  y(y+l)JI/i  / 


2 

(  «t  k\  ,  v-iA,l 

y+i 

Ll  +  *«. 

[l+ke,  */+yAfi*J 

(3.6) 

(3.7) 


P., 


2 _ y-1 

y+l  y(y+l)A/i 


(3.8) 


The  unknown  parameter  e2  appears  in  equations  (3.4)  and 

(3.7).  Therefore,  the  boundary  value  problem  of  the  differential 

equation  is  not  unique.  e2  is  dependent  upon  the  downstream  flow 

field.  In  reality,  it  cannot  be  predicted  and  will  require 

repeated  iterations  to  be  accurately  determined^^ .  Our 

objective  is  to  independentally  solve  the  stagnation  point 

equation.  In  practice  it  was  found  that  it  is  successful  to  use 

the  inviscid  shock  wave  to  approximate  e2 .  It  is  not  too 

different  from  the  stagnation  point  solution  obtained  by  using  a 

[12] 

linear  method  .  For  convenience,  a  small  amount  of  data  can 
be  used  to  prepare  an  e2~Ma9  curve  (See  Figure  2). 

In  equilibrium,  because  the  temperature  T  is  hidden  in  the 
equations,  the  shock  wave  equation  must  be  determined  by 
iteration . 
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Figure  1.  The  Coordinate  System 

1 .  shock  wave 

2 .  body 


Figure  2.  Curve 

IV.  Numerical  Method 

Because  some  physical  quantities  have  larger  gradients  near 
the  body  surface,  after  considering  the  stability  of  the 
numerical  method  and  the  accuracy  of  the  results,  we  used  a 
certain  coordinate  transformation  to  automatically  tighten  the 
lattice  at  places  where  the  gradient  is  large.  In  this  work,  we 
used  the  following  equation 


m  •*.  «*  ■*  »“  . '  ’  ’  ”  »  '  •  .  > 


.vV-vH 
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f=ln( 1+Gy) /InC 1+G)  (4.1) 

where  G  is  the  lattice  control  parameter.  When  the  Reynolds  / 

number  Re^  is  large,  G  is  chosen  to  be  larger.  When  Reo  is 
small,  G  is  chosen  to  be  smaller  depending  on  the  accuracy 
required . 

We  used  a  finite  difference  method  to  solve  the  boundary 
value  problem  of  the  normal  differential  equation  set.  Center 
difference  was  used  for  internal  points  to  ensure  that  the 
interrupt  error  is  of  the  second  order.  The  second  order 
equation  is  solved  by  a  chasing  method.  The  continuity  equation 
is  integrated  from  the  body  surface  to  the  shock  wave.  The 
equations  for  p0  and  p2  are  integrated  from  the  shock  wave  toward 
the  body  surface.  The  shock  wave  distance  can  be  determined  by 
two  methods:  one  is  to  the  Newtonian  Method  to  solve  non-linear 
equations  and  the  other  is  to  use  the  integral  of  the  continuity 
equation.  Numerical  calculations  show  that  the  latter  is  better 


than  the  former.  The  latter  can  improve  the  rate  of  convergence. 
The  formula  to  determine  the  shock  wave  using  the  Newtonian 


method  is 


—  a> 


PtV » 
d  ip-.V-J) 

de  *’ 


(4.2) 


where  ui  is  the  relaxation  factor  (O^io^l).  The  formula  based  on 
integrating  the  continuity  equation  is 


where 


F=2p0  /(Ui  -Kv0  )/(1+f  ’ ) 


(A) 


In  gas  equilibrium,  we  used  a  triple  sampling  function  to 
approximate  the  thermodynamic  function  table^ ^ ^ .  We  divided 
p  from  100*«>/1000  atm  into  one  section  and  other  pressure  as 
another  section.  We  used  8x17  and  22x37  to  approximate  the  table 
value.  The  maximum  relative  error  is  less  than  17.. 

V.  Results  and  Analysis 

We  calculated  over  a  wide  range  of  Mm  and  Reo  and  compared 
the  results  to  those  in  the  literature.  In  Figure  3,  we  compared 
the  calculated  ut  and  v0  with  those  obtained  using  a  linear 

[ii] 

method  .  The  computation  was  carried  out  for  Ma>=20,  Rea>=173, 
1730,  17300,  173000,  and  1730000;  y=1.4;  TW=300°K;  Pr=0.72;  u=K 
'/F~ ( 1+S/T)  . 


Figure  3.  The  u0  ,  vt  Cross-section 
1.  this  work 
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It  is  obvious  that  the  stagnation  point  solution  obtained  by 
using  the  e2  value  of  an  inviscid  flow  to  approximate  that  of  a 
viscous  flow  is  very  similar  to  the  solution  obtained  by  using  a 
linear  method.  For  a  quantitative  comparison,  Table  1  lists  the 
results  of  this  work  and  the  data  reported  in  reference  [13]. 

From  Table  1  we  can  see  that  the  maximum  error  occurs  in  ur  . 
However,  the  relative  error  is  less  than  0.5%.  If  the  effect  of 
downstream  flow  on  the  stagnation  region  is  not  considered  (i.e.e, 
=0),  then  the  maximum  error  in  the  shock  layer  may  reach  207., . 
Figure  4  shows  the  variation  of  ur  with  e2=0  and  e2^0. 

Figures  5  and  6  are  the  results  obtained  in  equilibrium. 

They  are  also  compared  to  the  viscous  equilibrium  gas  results 
obtained  using  a  linear  method  and  that  of  an  ideal  gas.  The 
agreement  is  very  good. 


Table  1  Stagnation  Data  Obtained  in  this  Work  and  by  Linear 

Method 


i  a  m 


S  *t  te 


vo  oo 


1. 

2. 

3. 

4. 

5. 


parameter 
method 
this  work 
linear  method 
this  work 
linear  method 
this  work 
linear  method 
this  work 
linear  method 


Figure  4.  Various 


Figure  5.  Gas  Equilibrium  Results 

1.  inviscid  flow 

2.  this  work 


Figure  6.  Gas  Equilibrium  Results 
1.  this  work 


1 .  this  work 

2.  inviscid  flow 


Figure  7  shows  a  comparison  of  body  surface  pressure 
obtained  by  using  this  method  to  that  calculated  from  an  ideal 
gas  stagnation  solution  over  a  wide  range  of  Mo.  Obviously, 
the  body  surface  pressure  of  a  viscous  flow  is  higher  than  that 
of  an  inviscid  flow  when  the  wall  is  cold.  This  means  that 
the  effect  of  viscosity  causes  the  body  surface  to  rise. 

Table  2  shows  the  variations  of  the  shock  wave  separation 
distance  c  and  body  surface  pressure  pw  with  Reynolds  number 
Rea  at  Mco=20.  With  increasing  Reynolds  number,  e  and  pw  are 
approaching  their  corresponding  inviscid  flow  values. 
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Table  2  Shock  Wave 


1 

P» 

MOO 

0.1274 

0.9214 

80000 

0.1237 

0.92057 

100000 

0.1298 

0.92055 

.  1000000 

0.13088 

'  0.92051 

oo 

0.1304 

0.32045 

1 .  parameter 

The  viscous  shock  layer  model  appears  to  be  applicable  over 
a  wide  range  of  Re^.  Based  on  our  analysis,  it  is  not  reasonable 
to  neglect  the  effect  of  downstream  flow  on  the  stagnation  point 
when  the  viscous  shock  layer  model  is  used  to  independently 
determine  the  solution  in  the  stagnation  region.  In  order  to 
overcome  this  shortcoming,  the  method  introduced  in  this  paper  is 
simple  and  effective. 
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NUMERICAL  METHOD  OF  VISCOUS  SHOCK-LAYER 
NEAR  THE  STAGNATION  LINE  OF  A  BLUNT  BODY 

Yu  Qingwen  Wang  Ruquan 
(The  Computing  Center,  Academia  Sini co) 

Abstract 

This  paper  gives  a  numerical  solution  of  the  viscous  shock-layer  equa¬ 
tions  described  the  flow  near  the  stagnation  line  of  a  blunt  body.  The 
difference  from  usual  treatments  is  that  authors  consider  the  influence  of 
down-stream  flow  on  the  stagnation  line  and  it  is  characterized  by  the 
boundary  conditions.  In  this  ease  an  unknown  parameter  which  depends 
on  the  shock  wave  is  approximated  by  inviscid  flow  data.  Numerical  re¬ 
sults  agree  with  those  given  by  the  solution  of  the  whole  flow  field  in 
front  of  the  blunt  body. 
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The  Perfection  and  Application  of  the  Flutter  Subcritical 

Response  Analysis 

Lu  Qizheng,  Lu  Guobao  and  Zeng  Weiqin 
(China  Aerodynamic  Research  and  Development  Center) 


I.  Introduction 

In  1979,  we  first  used  a  random  decrement  method  to  analyze 

[  1  ] 

the  subcritical  flutter  response  and  solved  the  fundamental 
problem  to  eliminate  the  effect  of  wind  tunnel  noise.  It  was 
pointed  out  in  reference  [1]  that  in  order  to  perfect  the 
analytical  method  and  to  put  it  into  practical  use  we  had  to:(1) 
reduce  running  time;  (2)  realize  on-line  analysis  and; (3) 
separate  modes  of  closely  spaced  frequencies.  To  this  end,  an 
on-line  analysis  system  was  built  between  1980  and  1981  to 
realize  on-line  analysis.  Furthermore,  the  running  time  of  a 
velocity  pressure  gradient  was  reduced  from  3  minutes  to  15 
seconds  and  a  curve  fitting  method  capable  of  separating  the 
subcritical  response  of  flutter  modes  of  two  closely  spaced 
frequencies  was  established.  In  addition,  the  accuracies  of 
measurement  and  analysis  of  the  subcritical  response  were 
investigated.  This  method  was  applied  to  various  models  with 
satisfactory  results. 


II.  On-Line  Analysis  System 

This  system  is  composed  of  a  tape  recorder,  a  power  filter, 
a  JCD-474  detector  and  a  DJS-21  computer.  After  each  run,  it 


performs  a  random  decrement  analysis.  By  single  mode  curve 
fitting,  the  frequency  f  and  relative  damping  coefficient  y  of 
the  flutter  are  given  to  instruct  the  continuation  of  the 
experiment.  In  addition,  a  power  spectrum  analysis  is  performed 
to  determine  whether  the  filter  band  selected  is  appropriate. 

The  so-called  random  decrement  method  involves  the  use  of  the 
same  voltage  U0  (representing  the  initial  displacement  of  the 
mode)  to  sample  N  sections  in  the  modal  response  time  process 
with  appropriate  filtering.  Then,  the  sum  is  averaged  to 
eliminate  the  effect  of  wind  tunnel  noise.  The  instantaneous 
response  of  the  mode  remains  to  subsequently  find  y  and  f. 

(1)  Sampling  Method 

In  1979,  when  the  analysis  was  done  on  the  Model  CF-700 
statistical  analysis  equipment,  because  we  did  not  use  a  special 
software  package,  samples  were  gathered  after  triggering  by  a 
constant  voltage  U0  .  Thus  a  larger  data  file  was  required.  The 
running  time  for  each  velocity  pressure  gradient  was  as  long  as  3 
minutes.  Currently,  we  are  using  an  overlapping  sampling  method 
which  treats  every  data  point  "greater  or  equal  to  U0 "  as  the 
beginning  of  a  sample.  In  addition,  the  raw  data  are  stored  in 
the  memory  of  the  computer  prior  to  performing  the  computation. 
Thus,  the  data  gathered  over  9  seconds  can  meet  the  accuracy 
requirement.  The  actual  running  time  is  reduced  to  15  seconds. 

(2)  Selection  of  Analytical  Parameters 

This  is  an  area  already  covered  in  reference  [1]  and  only 
some  additional  comments  are  made  here.  If  U„  is  too  small,  then 
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there  are 


This  paper  was  received  on  December  30,  1982  and  revised 
manuscript  was  received  on  August  27,  1983. 

more  intersects  in  the  response  time  period.  The  cumulative  / 6 

number  of  average  is  larger  which  favors  the  elimination  of  the 
effect  of  noise.  However,  because  the  average  amplitude  of 
response  is  lowered,  the  signal  to  noise  ratio  is  also  reduced. 
Therefore,  there  is  an  appropriate  range  of  U0  .  After  many 
trials,  U0  was  chosen  to  be  1.25  times  the  mean  square  root  of 
response  o.  In  addition,  comparing  using  "greater  or  equal  to  U0 
"  as  the  initial  condition  to  using  "equal  to  U0 "  as  the  initial 
condition,  we  found  that  N  could  be  increased  by  more  than  one 
fold  at  the  same  record  length.  This  also  favors  the  elimination 
of  the  effect  of  noise.  When  determining  the  number  of  samplings 
n  for  each  specimen,  in  addition  to  the  fact  that  the  random 
decrement  label  should  at  least  contain  6  waveform  cycles,  we 
must  also  consider  the  frequency  resolution  Af  requirement  for 
spectral  analysis  (Af=(1/nAt),  where  At  is  the  interval  between 
samplings.  The  larger  the  cumulative  average  number  N  is,  the 
smaller  the  effect  of  noise  becomes.  However,  the  length  of  the 
record  required  and  the  running  time  are  also  longer.  Therefore, 
it  is  necessary  to  choose  it  properly.  After  the  sampling  rate 
1/At  and  U0  are  selected,  the  length  of  the  sampling  time 
reflects  the  magnitude  of  the  cumulative  average  number.  It  was 
found  that,  when  U0=1.25o  and  (1/At)=1000  cycles/sec,  relatively 
stable  results  could  be  obtained  with  a  9  second  sampling  time. 
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(3)  Fitting  of  markers 

Characteristically,  the  initial  slope  of  every  marker  is 
zero.  Thus,  in  order  to  simplify  the  operation,  we  used  Ae~r"  cos  at 
for  curve  fitting.  When  the  signal  to  be  analyzed  only  contains 
a  single  mode,  such  a  fitting  can  result  in  relatively  accurate 
results.  If  the  signal  contains  two  modes,  although  this 
fitting  is  not  too  accurate,  yet  it  is  still  capable  of  reflect¬ 
ing  the  modal  variation  and  the  extent  of  approaching  the 
critical  flutter  point  with  varying  velocity  pressure.  It  is 
still  useful  in  experiments  for  online  monitoring.  As  for  more 
accurate  results,  a  power  function  curve  fitting  method 
suitable  for  two  modes  can  be  used  afterwards. 


III.  Curve  Fitting  Method 

In  order  to  process  the  random  decay  label  of  two  mixed 


modes  (also  possible  for  processing  single  mode),  a  curve  fitting 
program  written  in  FORTRAN  was  compiled  for  a  DJS-8  computer. 

The  important  features  of  this  method  are: 

A  complex  function  f  is  used  to  fit  a  label  y  from  the  DJS- 


21  computer. 


/“2 j  *'*'*''  (fl,  sin  a,t+b,  cos  o>,t) 


t=t^=KAt  (K=1,  2, . ,n) 

where  u>j  and  y .  are  the  characteristic  frequency  and  damping 
coefficient  of  the  mode  to  be  determined,  a^  and  b^  are  the 
coefficients  to  be  determined,  At  is  the  interval  between 
sampling,  and  K  is  the  sequence  of  the  data  points  in  the  label. 
For  a  dual  mode  system,  J=2.  The  set  of  parameters  to  be  solved 
is  (a*),  i— 1  , 2 . .  Let  3^  ~a^  ,  3]  ^  ,  a ,  ~  ,  a^  ~~  b  ^  ,  3,  —  ctj 


The  fitting  is  done  by  using  the  least  square  method  to 
minimize  the  residual  square  sum  of  the  fitting  function  f  and 


the  label  y 


J?-  2  Cy(ACA/)-/(A'A/)]» 


1  .e . , 


aR/aax  =  0  ( i=1 ,  2,  . ,m)  (2) 

Let  us  use  a  set  of  approximate  parameters  a^^and 


ai=ai  +Ai  (i=1>2> 


.WV  S*.  -.‘l  v”  * 


Then,  the  problem  to  find  la^}  becomes  the  problem  to  determine 
the  correction  terms  { Ai } .  By  expanding  f  around  a^^  into  a 
Taylor  series  and  neglecting  higher  terms  we  get 

'  fUK,  alt  a„ . A,  +  %i-  A,  + . +  -^-  A.  (4) 


where  **• 


““'“'"•■daT  r  da,  1  ' 

fK.=f(tK,  «»•*.  .«:•’> 

3/k»  ai»  gj.  . »g«> 

da,  da,  a.  =  a 


The  normal  equations  are  obtained  from  equation  (2). 

Furthermore,  we  introduced  a  damping  factor  v  into  diagonal 
elements  of  the  coefficient  matrix  in  order  to  relax  the 
selection  requirement  of  the  initial  values  a^®^  and  to  improve 
the  convergence  of  iteration.  The  equations  are  written  as 


where 


(a,,  +  »)  Ai  +  o,,  A,+  • 
a„  Ai+  (<*»«+») A,  +  • 


•  +  0j»  A«*®»»_ 

•  +  a,n 


(U,  Ai  +  o,j  A,+ . .  . + 

a  =  Y1  iZliL  .  iZ*S_f  (», 
a"  &  da,  da,  ’  ’ 

rr  o<*, 


(«,  /=  1, 2,  ***^ 


After  the  label  is  given,  's  are  first  selected.  Then, 

the  coefficient  matrix  and  the  right  hand  side  of  equation  (6) 
are  used  to  solve  Ai  by  using  the  Gaussian  elimination  method  in 

order  to  obtain  a^.  Let  £^  =  (A^/a^),  ( i= 1,2, . m).  The 

corrected  a^  is  used  as  the  initial  value  a^^  in  the  next 
iteration.  This  process  is  repeated  until  le^^e  (which  is  a 
pre-determined  small  number).  The  last  a^  obtained  is  the  system 
parameter . 

In  order  to  further  improve  the  converging  rate,  an  added 

optimal  step  factor  is  also  used.  If  the  1 iteration  does 

not  converge,  it  is  re-defined  as 

a  (l+1)_a  (1).A  (1)  n 

a^  =a^  +A^  .  Py  (.B; 

In  order  to  improve  the  accuracy,  the  label  was  standardized 
prior  to  fitting. 


IV.  Accuracies  of  Analysis  and  Measurement 
An  experimental  study  on  the  accuracy  of  the  subcritical 
flutter  response  was  carried  out  using  the  flutter  model  of  the  C2 
wing.  First,  we  investigated  the  statistical  accuracy  of  the 
subcritical  response  analysis.  The  signal  recorded  for  each  run 
was  analyzed  10  times  for  subcritical  response  and  the  mean  value  y 
and  mean  square  deviation  of  the  damping  coefficient  y  were 
obtained : 

v-4- i>. 


V  (V--  V>’ 


e 


(see  reference  [2]) 


In  addition,  if  we  assume  the  confidence  rate  a=0.95,  the 
confidence  region  for  y  is  defined  as  (7-t,  y+O: 

I— 

n 

‘  W3r  £<vT-l>7 


where  t.  can  be  looked  up  in  Appendix  4  in  reference  [2]  based  on 
the  degrees  of  freedom  K=n-1  and  the  confidence  level  a. 

Results  of  the  statistical  analysis  are  shown  in  Tables  1 
and  2.  The  y  and  f  for  model  Y 1-0000  were  obtained  by  an  on-line 
analysis  system.  The  y  and  f  for  model  Y2-DDDD  were  obtained 
from  curve  fitting  of  labels  recorded  by  an  on-line  analysis 
system.  In  both  examples,  the  standard  deviation  ay  of  the 
damping  coefficient  and  the  amplitude  of  the  confidence  region  t 
are  relatively  small,  except  for  the  second  order  values  of  the 
model  Y1-0000.  This  indicates  that  the  subcritical  response 
analysis  system,  including  on-line  analysis  and  curve  fitting, 
has  a  satisfactory  accuracy  to  meet  the  needs  of  model  testing. 


Table  1.  Statistical  Analysis  of  Accuracy  of  the  C2  Wing  Model 

(Y1-0000,  637  runs) 


M 

0  tkg/m1) 

/  ‘Hz,  | 

V  1 

Of 

t 

0.911  ! 

5204  ! 

66.19  I 

0.0658 

0.00334 

0.00244 

0.911 

5699  [ 

70.97 

0.0495 

0.00678 

0.00495 

0.911 

6197  ! 

72.78 

0.0237 

0.00183 

0.00133 

0.911 

6872 

73.69 

0.0164 

0.00120 

0.00120 

0.912 

6719  | 

74.48  1 

0.0061 

0.00036 

0  00027 

Table  2.  C2  Wing  Model  with  Missiles  (Y2-DDDD,  526  runs) 


M 

4  (kg/mM  | 

4 

j  (Hz) 

V 

Or 

t 

0.877 

3734 

■a 

0.01667 

0.000326 

0.000238 

0.881 

3913 

0.009S1 

0.000348 

0.000254 

0.888 

4069 

42.67 

0.00672 

0.000324 

0.000244 

0.894 

4809 

42.18 

0.00475 

0.000261 

0.000191 

0.898 

4493 

41.86 

0.00353 

0.000317 

0.000161 

0.982 

4683 

41.78 

0.00283 

0.000384 

0.000354 

For  the  same  model  state  (Y1-0000),  six  runs  were  performed 
in  repetition.  The  average  of  10  analyses  is  used  as  the 
"measured  value"  for  each  run.  The  least  square  method  is  used 
to  determine  the  critical  flutter  velocity  pressure  ,  by 
extrapolation.  The  relative  value  of  the  standard  deviation  of 
the  critical  velocity  pressure  between  runs  to  the  mean 

critical  velocity  pressure  q^  is  determined  as  follows: 

°V^P  =  1-30%  (c) 

The  relative  value  of  the  maximum  deviation  of  the  critical  /  70 

velocity  pressure  is 

qkpmax-?kp/qkp  =2-301  (A) 

which  indicates  that  the  reproducibility  of  the  measurement  can 
satisfy  the  requirement  of  the  model  experiment. 

V.  Application  Results 

(1)  Transonic  Flutter  Model  of  Model  A  Wing.  The 
experiment  was  carried  out  in  a  0.6  meter  transonic  wind  tunnel. 
Figure  1  shows  the  'r~q  curve  obtained  by  various  analysis  systems 
at  M=0.75.  An  after  the  fact  analysis  method  was  used  in  this 
experiment.  The  critical  velocity  pressure  obtained  by 
extrapolation  using  the  subcritical  method  differs  from  the 
directly  measured  value  by  less  than  8%. 
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Figure  1.  Relative  Damping  Coefficient  vs.  Velocity  Pressure 

for  a  Model  A  Transonic  Wing  (M=0.75) 

1.  oCF-700  subcritical  measured  points 

2.  aDJS-21  subcritical  measured  points 

3.  Vcurve  fitting 

4.  (DJS-8)  subcritical  measured  points 

5.  directly  measured  critical  points 

(2)  Low  Velocity  Flutter  for  Model  B  Aircraft.  The 
experiment  was  carried  out  in  a  4mx3m  low  velocity  wind  tunnel. 

In  the  case  of  a  15  ton  commercial  load,  0303  hoist,  no  fuel,  and 
with  symmetric  flutter,  because  the  critical  flutter  point  does 
not  appear  up  to  the  allowed  wind  speed  of  58m/sec,  we  can  only 
rely  on  extrapolation  based  on  the  subcritical  method.  For  both 
modal  states,  all  labels  showed  a  beat  resonance  effect  due  to 
the  mixing  of  two  closely  spaced  frequency  modes.  Therefore,  it 
is  necessary  to  go  through  curve  fitting.  The  resultant  -y-q 
curve  for  the  flutter  mode  is  shown  in  Figure  2.  We  can  see  that 
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the  critical  flutter  velocity  determined  with  the  subcritical 
method  is  in  good  agreement  with  the  directly  measured  value  and 
the  theoretical  value. 

(3)  Supersonic  Flutter  Model  of  60°  Delta  Flat  Wing.  The 
experiment  was  performed  in  a  0.6  meter  supersonic  wind  tunnel. 
The  M  numbers  were  1.5,  2.0  and  2.5.'  It  was  discovered  in  the 
experiment  that  the  model  exhibited  a  near  sinusoidal,  but 
slightly  damped  vibration  over  a  wide  velocity  pressure  range. 

It  was  difficult  to  directly  measure  the  critical  point.  For 
this  reason,  the  subcritical  response  y~c\  curve  shown  in  Figure  3 
was  used  to  determine  the  critical  point  by  defining  the  turning 
point  at  which  y  becomes  small  (-*<0.01)  and  flat. 


Figure  2.  Relative  Damping  Coefficient  vs.  Velocity  for  Low 
Velocity  Model  B  Wing 

1.  subcritical  measured  point 

2.  directly  measured  critical  point 

3.  subcritical  measured  point 

4.  directly  measured  critical  point 


61 


(4)  Transonic  Flutter  Model  C  Wing.  The  experiment  was 
carried  out  in  a  0.6  meter  transonic  wind  tunnel.  Typical 
results  are  shown  in  Figure  4.  In  this  experiment,  an  on-line 
analysis  on  the  subcritical  response  was  performed.  The  critical 
points  obtained  by  extrapolation  using  the  subcritical  method  are 
within  27,  of  directly  measured  values.  As  compared  to  the  87, 
deviation  in  the  off-line  analysis  of  Model  A,  because  on-line 
analysis  enables  the  measured  point  to  get  closer  to  the  critical 
point,  extrapolated  results  are  more  accurate. 


4000  5000  0000  7000  1000  1 

Figure  3.  Relative  Damping  Coefficient  vs.  Velocity  Pressure  for 
a  Delta  Flat  Supersonic  Wing 

1.  determined  critical  point 


Figure  4.  Relative  Damping  Coefficient  vs.  Velocity  Pressure 
for  Model  Transonic  Wing  M  =  0.866 

1.  subcritical  measured  point 

2.  directly  measured  critical  point 

(5)  Transonic  Flutter  of  C£  Wing  with  Missiles.  The 
experiment  was  carried  out  in  a  0.6  meter  transonic  wind  tunnel. 
Because  the  subcritical  method  was  used,  more  experimental 
results  were  obtained  with  each  model,  including  the  comparison 
of  the  critical  flutter  velocity  pressure  at  various  hanging 
positions  and  the  effect  of  transonic  compressibility  of  Y2-DDDD. 
This  fully  demonstrates  the  superiority  of  this  method.  The 
subcritical  analysis  of  the  typical  Y2-DDDD  model  is  shown  in 
Figure  5. 
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Figure  5.  Relative  Damping  Coefficient  vs.  Velocity  Pressure 
for  Model  C2  Wing  with  Missiles  (M=0.866) 

1.  single  mode  fitting 

2.  dual  mode  fitting 

(a)  y-q  curve 

(b)  fUq  curve 


VI.  Conclusions 

The  subcritical  flutter  response  analysis  based  primarily  on 
a  random  decrement  with  power  spectral  analysis  is  a  successful 
one.  The  power  spectrum  method  is  effective  in  determining  the 
mode  of  analysis  and  filter  bandwidth.  In  the  event  that  two 
closely  spaced  frequencies  are  mixed,  curve  fitting  is 
imperative.  The  establishment  of  an  on-line  analysis  system  can 
improve  the  accuracy  of  measurement  and  reduce  the  running  time. 
This  method  can  be  used  in  transonic,  supersonic  and  low  speed 
flutter  experiments. 
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THE  PERFECTION  AND  APPLICATION  OF  THE  FLUTTER 
SUBCRITICAL  RESPONSE  ANALYTICAL  METHOD 

Lu  Qizheng  Lu  Guobao  Zeng  Weiqin 
(Chine  Aerodynamic  Research  ft  Development  Center) 

.  ..  Abstract  *—  -  ~ 

The  Random  Decrement/Power  Spectrum  subcritical  response  analytic*! 
method  was  used  recently  in  the  wind  tunnel  model  flutter  test.  This  pa¬ 
per  describes  the  perfection  and  application  of  the  method,  and  includes 
the  establishment  of  the  on-line  analytical  system,  the  establishment  of 
curve  fitting  method  with  exponential  functions,  the  study  of  accuracy, 
and  the  results  from  application  of  the  method.  The  method  has  foilwing 
features,  the  effect  of  the  wind  tunnel  flow  noise  can  be  eliminated  ef¬ 
fectively,  the  running  time  is  shorter,  the  modes  of  closely  space  freq¬ 
uencies  can  be  identified.  Its  feasibility  is  good.  The  practical  application 
in  the  tests  of  Model  A,  B,  C  and  a  60°  delta  flat  wing  model  shows  that 
the  method  can  be  applied  to  transonic,  supersonic  and  low  speed  wind 
tunnel  flutter  tests. 


Calculation  of  Octagonal  Wall  Interference  Factor  Using 

Conformal  Mapping 

(Northwestern  Polytechnical  University) 

I.  Introduction 

Currently,  the  cross-section  of  many  low  velocity  wind 
tunnels  is  octagonal.  Furthermore,  the  6  curve  given  in 
references  [3]  and  [6]  is  used  for  wall  correction. 

There  are  primarily  two  methods  to  determine  the  value  of  6 
one  is  the  image  method  which  is  capable  of  calculating  the  6 
values  for  circles,  squares,  rectangles,  ellipses  and  the  other 
is  conformal  mapping  which  is  capable  of  calculating  the  values 
of  6  for  rectangles,  squares,  ellipses  and  octagons.  However, 
based  on  information  available  to  date,  we  still  could  not  find 
an  accurate  formula  to  calculate  the  interference  factor  for  an 
octagonal  wall  based  on  the  two  methods  mentioned  above.  The 
methods  used  in  papers  discussing  this  problem  are  appropriate. 
The  method  in  this  work,  however,  can  provide  an  accurate  form”! 
for  calculation.  In  addition,  this  method  is  also  much  simpler 
than  those  approximation  methods.  It  is  capable  of  calculating 
the  6  values  for  rectangles,  squares  and  octagonal  wall 
interference  factors. 

In  the  following  calculation,  let  us  assume  that  the  wing 
load  is  uniform,  i.e.,  the  wing  is  replaced  by  a  pair  of 
vortices.  In  addition,  it  also  assumes  that  the  wing  is  located 
in  the  plane  of  symmetry  in  the  test  section  of  the  wind  tunnel. 


II.  General  Equation  for  Calculating  Wall  Interference 

Factor  of  Any  Shape  Using  Conformal  Mapping.  The 

interference  of  a  circular-section  wall  can  be  easily  resolved. 

It  is  natural  to  think  that  whether  this  problem  can  be  resolved 

if  the  conformal  mapping  formula  for  transforming  the  inside  of 

an  arbitrary  wind  tunnel  to  a  circular-section  is  known. 

This  problem  had  been  proven  in  reference  [1] . 

Let  us  assume  that  the  following  transformation  formula 

€  =  f(Z)  (2.1) 

can  transform  the  cross-section  of  a  specific  wind  tunnel  into  a 
unit  circle  |Z|-1  on  the  Z  plane  by  conformal  mapping. 
Furthermore,  the  mapping  is  conformal  for  every  point.  Let  us 
assume  that  the  wing  in  the  tunnel  is  replaced  by  a  pair  of 
vortices  (r,-r)  spaced  by  b  (the  wing  span),  as  shown  in  Figure 
1. 

If  the  vortices  are  located  at  %x  and  £2  on  the  £  plane, 

they  are  transformed  to  Zx  and  Za  on  the  Z  plane  without  changing 

M] 

its  intensity 

Based  on  reference  Cl],  when  the  wing  is  located  on  the  real 
axis  and  the  center  of  the  wing  span  is  at  the  origin,  if  the 
interference  angle  at  the  center  of  the  wing  span  is  used  to 
represent  the  value  on  the  entire  wing,  then  the  value  of  6  is: 

—  *  (J— _LU  (2.2) 

2 sb*  \  2/  (0)  \  l  2/1 

This  paper  was  received  on  August  18,  1983  and  revised  on  October 
5. 
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where  Ay  is  the  cross-sectional  area  of  the  experimental  section 
of  the  wind  tunnel. 

When  the  span  is  large  relative  to  the  dimension  of  the 
experimental  section  of  the  wind  tunnel,  it  is  more  appropriate 
to  use  the  average  interference  angle  along  the  wing  span.  In 
this  case,  we  get: 


Aw 


~i  log 


,,  (  l  \ 

(2  l\l 

f  \  t) 

\  l  2/2 

(2.3) 


III.  Conformal  Mapping  Formula  from  Unit  Circle  to  Octagon 
It  was  given  in  reference  [2]  that 


where  a^(K=l,2, . ,n)  -  apex  angles  of  the  polygon  in  the  unit 

of  ti . 

aK(K=1,2, . ,n)  -  points  on  the  unit  circle 


corresponding  to  apices. 

From  Figure  2  one  can  see  that  A',  B',  C',  D',  E',  F',  G' 
and  H'  on  the  octagon  correspond  to  points  A,  B,  C,  D,  E,  F,  G 
and  H  on  the  unit  circle. 

P°int  A  B  C  ■>-  D  E  F  -  G  H _ 

a„  «• e‘‘#*  «*i#* 


2. 


Figure  2.  Relation  Between  5  Plane  and  Z  Plane 

1.  €  plane 

2.  Z  plane 


From  Figure  2  one  can  see  that  the  value  of  is  as 
follows : 


Z.A'—a*,  o+a'»—-, 

z 


By  substituting  the  and  obtained  into  equation  (3.1) 
we  get 


(3.2) 


Equation  (3.2)  includes  three  parameters:  0X  ,  0a  and  a. 

For  a  specific  wind  tunnel,  these  three  parameters  are  uniquely 
defined. 

When  a  is  fixed,  two  length  ratios  (0 ' Q ' /0 1 P ' )=A  and 
(A'P'/O'P' )=B,  can  be  used  to  define  an  arbitrary  octagonal 
section. 

1.  Calculation  Formula  for  O'P' 

From  Figure  2  one  can  see  that  in  this  case,  Z=x.  Because 
it  is  a  unit  circle,  equation  (3.2)  becomes: 

C  <1 +  ideas’  0J*-,[(l  +  *O,-4*'cos’0, 

Let  a=  3/4  +  1/4B,  B=[1,-1).  Then 
ct-1=-1-B/4 ,  1 / 2-a=-1+B / 4 

Substituting  it  into  the  above  equation,  we  get 

=  +  — 4jf’  cos  51]*1'T£C(l  +  Jc,),-4al  cos*  0,] '^dx 

2.  Calculation  Formula  for  O'Q' 

From  Figure  2  one  can  see  that  Z=iy.  For  convenience,  let 

y=x.  After  substituting  it  into  equation  (3.2),  we  get 

cos'  0,]‘iT,Cn-x,),  +  4jc,cos,flJ]'iT1dx  (3.4) 

3.  Calculation  Formula  for  A ' P  * 

i  9 

From  Figure  2  one  can  see  that  Z=e  .  After  substituting  it 
into  equation  (3.2)  we  get 

'  Cs«n  (0,  —  0)sin(0  +  0i>3"-^Tsin(0,—  0)sin(0,  +  0)3  dB 


(C) 

(3.3) 


(D) 


On  the  above  integral  is  a  generalized  one  when  6-*©!  .  It 
should  be  treated  properly.  For  this  purpose,  we  introduced  a 
new  variable,  i.e., 

( Qt  -0 ) £ , then  d£=3+B/4  (ex -e)-1-B/4de 
Therefore , 

CL  («.-*,  J-4*.  -iy  ]'^  dl 

After  substituting  it  into  the  expression  for  A'P',  we  get: 


(3.5) 


(E) 

(F) 


IV.  Special  Cases 

1 .  Formula  for  Normal  Octagonal  Section 

The  formula  for  a  normal  octagonal  section  can  be  obtained 
by  letting  0^22.5°  and  0a=67.5°  in  the  above  equation. 

2.  Formula  for  Rectangular  Section 

The  formula  for  rectangular  section  can  be  obtained  by 
letting  0X  =0a:  in  the  above  equation.  In  addition,  only  the 
integrals  of  O'P'  and  0 ’ Q '  are  required  to  be  calculated  when 
calculating  the  interference  factor  of  a  rectangular  section 
because  A ' P ' =0 ' Q ' . 


(1)  Formula  for  O'P' 


In  this  case,  B=- 1 , 0,.  =0a  =0 .  By  substituting  these  factors 

f  f>  .  ,dx 

into  equation  (3.3)  we  get 
(2)  formula  for  O'Q' 


4.1 


In  this  case  0^=o0 =0  and  P  =  -1 .  Substituting  these  into 
equation  (3.4)  we  get 

dx 


OQ'  =  \ 


»\/  <1  —  4X;  COSJ 


(4.2) 


/  8 1 
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3.  Formula  for  square  section 

When  0'P'=0'Q',  i.e.,  a  square,  it  is  determined  by  the 


following  formula: 


dx 

V  i-i-jcr 


(4.3) 


V.  Conclusions 

1.  When  K<0.6,  the  6  value  of  an  octagonal  section 
calculated  by  using  this  method  is  very  close  to  that  of  an 
elliptical  section  of  the  same  height  to  width  ratio.  When  K> 
0.6,  the  difference  is  larger,  as  shown  in  Figure  3. 

2.  The  6  value  of  a  regular  octagon  calculated  by  using 
this  method  is  identical  to  those  reported  in  the  literature,  as 
shown  in  Figure  4. 


Figure  3.  Comparison  of  6  Values  of  Octagonal  Sections 

(x=wind  tunnel  height/width  =0.778  and  k=wing  span/ 
wind  tunnel  width) 

1.  equation  (2.3) 

2.  ellipse  in  reference  [4] 

3.  equation  (2.2) 
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Figure  4.  Comparison  of  6  Values  of  Octagonal  Sections 


1.  circle  in  reference  [4] 

2.  equation  (2.2) 

3.  equation  (2.3) 

3.  It  is  suggested  that  this  method  be  used  to  calculate 
the  corresponding  curve  for  a  given  octagon,  instead  of  the  curve 
derived  in  references  [3]  and  C6D. 
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WALL-INTERFERENCE  CALCULATION  OF  WIND 
TUNNEL  WITH  OCTAGONAL  SECTIONS  USING 
CONFORMAL  MAPPING  METHOD 

Xia  Yushun  and  Lin  Chaoqiang 
(.Northwestern  Polytechnical  University) 

Abstract 

.  The  conformal  mapping  formula  is  used  for  the  wall-interference  calcul¬ 
ation  of  wind  tunnel  with  octagonal  sections.  The  parameters  in  the 
mapping  formula  can  be  easily  determined  by  computer.  As  particular 
examples,  the  results  for  rectangular,  square  and  regularoetagon  sections 
are  also  given  in  closed  form.  Some  typical  results  are  plotted  and  com¬ 
pared  with  other  results. 
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Cryogenic  High  Reynolds  Number  Transonic  Wind  Tunnel  with 
Pre-cooled  and  Restricted  Flow 


/  87 


Pan  Ruikang 

(China  Aerodynamic  Research  and  Development  Center) 

Introduction 

Since  the  sixties,  many  foreign  countries  were  concerned 
about  the  problem  of  insufficiently  low  Reynolds  numbers  in  wind 
tunnel  experiments  which  seriously  affected  the  performance  of 
aircraft.  Due  to  insufficiently  low  Reynolds  numbers,  areas  in 
aerodynamics  primarily  caused  by  viscosity  and  certain 
aerodynamic  phenomena  with  strong  coupling  of  viscous  surface 
flow  and  inviscid  external  flow  cannot  be  simulated.  Therefore, 
aerodynamic  researchers  proposed  that  the  entire  range  of 
Reynolds  number  encountered  in  flight  should  be  simulated  in  the 
wind  tunnel.  NATO  countries  and  the  U.S.  organized  special  task 
teams  to  study  this  problem.  Several  high  Reynolds  number 
transonic  wind  tunnels  have  been  developed,  including  the  low 
temperature  nitrogen  wind  tunnel,  Ludvic  tube  wind  tunnel,  jet 
wind  tunnel  and  Evans  wind  tunnel. 

Historically,  Osborne  Reynals  began  by  focusing  his 
attention  on  the  effect  of  low  temperature  gas  on  flow 
characteristics.  In  addition,  low  temperature  experiments  were 
done  and  a  suggestion  to  lower  air  temperature  was  made.  In 
1920,  William  Margoulis  (NACA)  was  the  first  one  to  introduce  the 
concept  of  low  temperature  wind  tunnels  to  consider  the  effect  of 
pressure,  temperature  and  the  experimental  gas  on  the  aerodynamic 
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experiments.  In  1945,  in  order  to  reduce  the  power  of  high  speed 

wind  tunnels  and  to  improve  the  economics,  Smelt  suggested  that 

the  gas  and  temperature  should  be  appropriately  chosen.  Due  to 

various  reasons,  low  temperature  gas  did  not  receive  any 

Cl] 

attention  .  In  1977,  Storllery  in  England  introduced  the  idea 

of  using  the  adiabatic  expansion  of  air  stored  in  a  high  pressure 

container  to  obtain  the  low  temperature  desired.  In  the 

experiment,  the  control  of  a  piston  movement  is  used  to  maintain 

[2] 

the  temperature  and  pressure  at  the  stagnation  point  .  In 
1979,  Curt  Nelander  of  Sweden  introduced  a  new  scheme  in  the 
first  international  discussion  meeting  of  low  temperature  wind 
tunnels.  In  his  scheme,  a  compressor  is  used  to  compress  the  gas 
exhausted  from  the  wind  tunnel  which  is  stored  in  the  low 
temperature  pipe  underground.  A  portion  of  the  compressed  air  is 
stored  in  a  high  pressure  vessel  and  another  portion, 
approximately  387o,  passes  through  an  air  cooling  device.  In  the 
experiment,  the  air  in  the  high  pressure  vessel,  after  adiabatic 
expansion  and  throttling  is  mixed  with  low  temperature  air  to 
ensure  that  the  required  temperature  and  pressure  are  attained  in 
the  stable  section  of  the  wind  tunnel.  In  1980,  Lian  Qixing  of 
Beijing  Institute  of  Aeronautics  and  Astronautics  introduced  a 
plan  to  use  heat  exchangers  for  cooling^^.  In  this  scheme,  the 
air  in  the  high  pressure  storage  system  is  adiabatically  expanded 
to  allow  the  heat  reservoir  to  attain  the  desired  low 
temperature.  In  this  experiment,  medium  pressurized  air  passed 
through  the  heat  reservoir  to  reach  the  desired  low  temperature . 
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Since  1975,  we  surveyed  various  high  Reynolds  number,  Re, 
wind  tunnels.  Primarily  we  believe  that  the  cryogenic  nitrogen 
wind  tunnel,  particularly  continuous  low  temperature  nitrogen 
wind  tunnel  such  as  NTF  of  NASA,  is  more  desirable  among  various 
high  Re  number  wind  tunnels.  Because  of  the  high  investment  of 
building  such  a  wind  tunnel,  it  is  significant  to  seek  new 
technical  approaches.  This  paper  describes  one  such  approach. 

Although  the  cryogenic  high  Re  number  wind  tunnel  with  pre¬ 
cooled  and  restricted  flow  introduced  in  this  paper  still  have 
some  technical  difficulties  in  high  pressure  low  temperature 
engineering,  however,  it  is  relatively  easier  as  compared  to 
other  cryogenic  air  wind  tunnels.  Because  of  the  pre-cooling  and 
flow  restriction,  not  only  the  experience  gathered  in  the 
construction  and  operation  of  the  gas  storage  systems  of 
hypersonic  wind  tunnels  and  the  experience  acquired  in  the  use  of 
intermittent  wind  tunnels  can  be  applied,  but  also  the  cost  of 
building  such  wind  tunnels  is  much  lower  than  that  of  a  nitrogen 
cryogenic  wind  tunnel.  Based  on  the  pre-cooling  and  flow 

Manuscript  received  on  August  18,  1982,  revised  manuscript 
received  on  October  26,  1983. 

restriction  approach,  this  paper  introduces  a  2.4m  high  Re  number/88 
cryogenic  transonic  wind  tunnel. 

Cooling  System  for  Pre-cooling  and  Flow  Restriction 

The  temperature  of  air  can  be  lowered  by  heat  transfer 
through  a  medium  or  by  refrigeration^^.  Because  a  pressure 
regulating  valve  is  used  to  exhaust  a  large  amount  of  stored  air 
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in  an  intermittent  wind  tunnel,  therefore,  the  experimental 
process  itself  is  a  refrigerator  process  for  air.  The  adiabatic 
expansion  of  stored  gas  and  restriction  of  flow  by  the  pressure 
regulating  valve  is  used  to  refrigerate  the  air  to  a  low 
temperature.  Although  the  temperature  lowering  effect  due  to 
flow  restriction  of  the  valve  can  be  calculated,  however,  it  is 
more  convenient  to  look  it  up  in  the  temperature-entropy  plot 
from  an  engineering  point  of  view^^’^^.  Because  the  pressure 
of  stored  gas  in  a  transonic  wind  tunnel  is  relatively  low,  the 
temperature  drop  due  to  the  integrated  flow  restriction  effect  is 
small.  However,  the  temperature  drop  caused  by  the  adiabatic 
expansion  in  the  gas  storage  vessel  is  large.  In  order  to  avoid 
this  temperature  drop  effect,  a  heat  reservoir  is  stored  in  the 
container  to  control  the  temperature  drop. 

The  adiabatic  expansion  or  air  in  a  container  or  flow 
restriction  by  a  throttle,  or  both  can  be  used  to  refrigerate  the 
cryogenic  wind  tunnel  itself.  Plans  introduced  by  Stollery  and 
Nelander  require  the  use  of  complicated  systems  to  control 
fluctuations  of  temperature  and  pressure  at  the  stagnation  point. 
In  this  aspect,  pre-cooling  with  flow  restriction  and  adiabatic 
expansion  are  relatively  simple.  The  following  is  a  comparison 
of  these  two  plans. 

Figure  1  shows  an  adiabatic  expansion  cooling  system.  Air 
from  the  atmosphere  or  in  10  atm  containers  is  compressed  to  175 
atm  and  stored  in  vessels.  The  valve  A  is  opened  to  allow  the 
air  from  the  high  pressure  vessel  to  expand  adiabatical ly  to 
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lower  the  pressure  from  175  atm  to  10  atm.  Then,  valve  A  is 
closed.  In  the  experiment,  the  flow  restricting  valve  controlled 
the  air  to  reach  state  d  (5  atm,  154  atm). 


Figure  1.  The  Adiabatic  Expansion  Air  Cooling  System  and  the 
Temperature  vs .  Entropy  Curve 


1. 

2. 

compressor 

drier 
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gas  storage 
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Figure  2.  The  Pre-cooled  and  Flow  Restricted  Air  Refrigeration 
System  and  the  Temperature  vs.  Entropy  Curve 

1.  cooling  unit 

2.  gas  storage  vessel 

3.  temperature  stabilizer 

4.  pressure  stabilizer 

5.  isobaric  throttle 

6.  pressure  regulating  valve 

7.  hypersonic  wind  tunnel 

8.  compression 

9.  cooling 

10.  flow  restriction 

Figure  2  shows  the  pre-cooled,  flow-restricted  refrigeration 


system,  the  gas  undergoes  the  processes  of  compression  Ca-b), 
cooling  (b-c),  isothermal  pressure  reduction  (c-c1 2 3 4 5 6 7 8 9 10)  and  flow 
restriction  (e'-d).  The  compressor  compresses  the  air  to  220 


atm.  It  is  then  cooled  to  215K  by  passing  through  the 
refrigerator.  It  is  finally  stored  in  a  high  pressure  vessel. 

In  the  experiment,  the  pressure  regulating  valve  and  the  isobaric 
flow  throttle  are  opened  simultaneously.  The  former  keeps  the 
pressure  stabilizing  box  at  150  atm  at  all  times  and  finally 
maintains  the  stagnation  pressure  in  the  stabilized  section  in 
state  d.  The  pre-cooling  system  in  the  system  is  a  conventional 
cooling  unit,  which  may  be  a  freon  refrigeration  system  or  a 
flow-restricted  pre-cooling  system  or  an  exhaust  gas  pre-cooling 
system. 

The  choice  of  either  refrigeration  plan  should  be  made  based 

on  technical  difficulties  and  experimental  requirements.  It 

should  also  take  factors  such  as  energy  consumption,  capital 

cost,  and  experimental  expenses  into  account.  Table  1  shows  a 

comparison  of  these  two  systems  in  a  transonic  wind  tunnel  with  a 

2.4m  x  2.4m  experimental  cross-section.  From  the  table  one  can 

see  that  the  storage  vessel  in  the  air  expansion  system  is  3  x  10 

3 

m  which  is  impossible  to  materialize  in  practice.  Although 
Stollery's  plan  was  improved,  however,  it  still  requires  a 
considerably  large  volume.  Thus,  it  is  difficult  to  control  the 
piston.  Nelander's  plan  was  indeed  a  significant  improvement. 
However,  it  is  difficult  to  control  the  temperature  and  pressure 
of  the  flow  in  the  stable  section.  On  the  contrary,  pre-cooled 
flow-restricted  refrigeration  systems,  no  matter  whether  pre¬ 
cooled  or  temperature  and  pressure  controlled,  can  be  realized 
from  an  engineering  point  of  view.  The  gas  supply  system  is 


similar  to  that  of  a  hypersonic  wind  tunnel.  In  this  case,  this 
system  can  be  shared  with  a  hypersonic  wind  tunnel.  Not  only  the 
technical  difficulty  is  low,  but  also  the  capital  investment  it 
less.  It  is  a  more  rational  plan. 


Table  1.  Comparison  of  Adiabatic  Expansion  and  Pre-cooled 
Slow-restricted  Systems 
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1.  item 

2.  adiabatic  expansion  refrigeration  system 

3.  pre-cooled  flow-restricted  refrigeration  system 

4.  gas  storage  condition 

5.  gas  pressure  atm 

6.  medium  pressure  atm 

7.  final  pressure  atm 

8.  initial  temperature  K 

9.  lowest  temperature  K 

10.  weight  of  stored  gas  T 

11.  volume  of  stored  gas  m* 

12.  residual  stored  gas  T 

13.  temperature  drop  in  experiment  K 

14.  temperature  compensation  K 

15.  weight  of  aluminum  heat  reservoir 

16.  energy  consumed  by  system  Kg.M 

17.  major  similar  and  different  equipment  and  systems 

18.  pre-cooling  system 

19.  drier 

20.  volume  of  storage  vessel  for  gas  exhausted  in 
lowering  temperature 

21.  pressure  stabilizing  box 

22.  high  speed  high  pressure  valve 

23.  isobar ic  throttle 

24.  to  ensure  the  humidity  in  air  is  0.002gH90/Kg  air 

25.  300,000  Z 

26.  none 

27.  yes 

28.  none 

29.  cooling  air  from  300K  to  215K 

30.  maintaining  150  atm  during  experiment 

31.  yes 

32.  none 

33.  yes 


The  2.4  Meter  High  Re  Number  Transonic  Cryogenic  Wind  Tunnel 
The  dimensions  of  the  experimental  section  of  the 
intermittent  high  Re  number  Cryogenic  Air  Wind  Tunnel  are  2.4m  x 
2.4m.  The  M  number  in  the  experiment  is  0.5^1. 6.  It  is  required 
that  within  the  range  of  M  encountered  in  the  experiment  that  the 
Re  number  is  4  x  10^  with  respect  to  a  0.24m  characteristic 
length.  That  means  Re=16.7  x  10^  with  respect  to  a 
characteristic  length  is  one  meter. 
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In  order  to  realize  this  experimental  state,  stagnation 
pressure  and  temperatures  at  various  M  numbers  are  chosen  from 
Figure  3.  If  the  design  state  is  M=0.8,  then  the  corresponding 
stagnation  pressure  P0 =5  atm  and  stagnation  temperature  T0=154K. 

A  wind  tunnel  is  composed  of  a  gas  supply  system,  the  main  / 90 
wind  tunnel  and  a  waste  gas  recovery  system.  Figure  4  shows  the 
variation  of  states  in  the  process  a-b-c-c'-d  in  the  wind  tunnel 
experiment.  In  reality,  c-c'  does  not  exist.  Instead,  it  varies 
according  to  c-c".  In  this  case,  the  final  temperature  may  reach 
138K.  From  an  engineering  point  of  view,  if  we  assume  that  the 
temperature  loss  along  the  path  is  equal  to  the  cc"  flow- 
restricted  temperature,  then  c-c'  is  selected.  The  process  c'-d 
is  realized  by  the  isobaric  flow  restricting  throttle.  The  final 
point  d  corresponds  to  the  gas  parameters.  P„  =5  atm  and  T0=154K 
in  the  experimental  section.  The  total  picture  of  the  wind 
tunnel  is  shown  in  Figure  5. 


Figure  3.  Stagnation  Point  and  Temperature  of  Cryogenic  Air 

Wind  Tunnel 

1.  simulating  Re  =16.7  x  10^ 

2.  avoid  air  condensation 

In  the  gas  supply  system,  there  are  five  Model  H22  axi-flow 
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air  compressors.  The  inlet  air  flow  is  130m  /mm  per  unit  and  the 
exhaust  pressure  is  220  atm.  The  axial  power  N=630KW.  The  air 
filling  time  for  the  first  experiment  is  5  hours.  Under  designed 
conditions,  the  filling  time  is  2  hours.  The  pre-cooling  system 
not  only  provides  cooled  air  but  also  makes  the  humidity  in  the 
air  to  be  under  0.2g  water /kg  air.  Two  sets  of  heat  exchangers 
are  used  alternately.  When  one  unit  is  operating,  the  other  unit 
is  being  defrosted.  The  volume  of  the  high  pressure  gas  storage 
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vessel  is  1290m  .  It  is  made  of  09Mn2VR  low  temperature  steel. 

The  temperature  stabilizing  device  is  used  to  compensate  the 
temperature  drop  caused  by  the  adiabatic  expansion  of  gas  in  the 
high  pressure  vessel.  The  initial  temperature  drop  is  estimated 
to  be  18. IK.  It  is  necessary  to  compensate  8. IK.  Thus,  it  is 
required  to  place  530T  of  LF5  aluminum  alloy  in  the  temperature 
stabilizer.  After  considering  the  actual  conduction  effect,  it 
should  be  increased  by  207o,  i.e.,  636T.  When  designing  the  heat 
exchanger  element,  the  contact  area  should  be  increased  and  the 
wall  thickness  should  be  decreased. 

The  pressure  regulating  valve  and  the  isobaric  throttle  are 
important  components  of  the  system.  The  former  can  maintain  the 
pressure  in  the  pressure  stabilizing  box  and  the  error  is 
controlled  within  the  range  of  0. 5-1.0%.  The  latter  should 
maintain  the  pressure  drop  at  a  given  M  number  in  the  experiment. 
Under  the  condition  that  the  flow  cross-section  of  the  valve 
remains  unchanged,  the  pressure  fluctuation  in  the  stable  section 
should  be  less  than  +0.57«.  The  range  of  adjustment  of  the 
pressure  ratio  of  the  pressure  regulating  valve  is  0.68^1.  The 
inlet  diameter  is  2m  and  the  outlet  diameter  is  2.4m.  The 
pressure  ratio  adjustment  range  of  the  isobaric  flow  restricting 
valve  is  0.06~0.01.  The  inlet  diameter  is  2.4m.  The  pressure 
regulating  ranges  of  these  two  valves  are  narrow.  It  is  easy  to 
realize  flow  restriction.  It  is  feasible  to  use  the  common  hole 
or  slot  acoustic  damper  valves. 

The  end  of  the  exhaust  diffusion  section  of  the  wind  tunnel 
is  equipped  with  a  wave  stabilizer  which  is  used  to  control  the 
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position  of  the  wave  to  allow  a  stable  M  number  in  the 
experimental  section.  In  the  experiment,  the  pressure  in  the  gas 
collecting  vessel  is  increasing.  The  shock  wave  in  the  wave 
stabilizer  continuously  moves  forward.  However,  the  wave 
stabilizer  can  stabilize  the  wave  downstream  from  the  throat  to 
maintain  the  aerodynamic  condition  for  the  experiment. 

The  exhaust  gas  collecting  vessel  collects  the  low 
temperature  exhaust  gas.  A  part  of  it  is  the  coolant  of  the  pre¬ 
cooling  system  and  the  other  part  is  sucked  into  the  compressor 
to  be  re-cycled.  For  this  reason,  the  pressure  of  the  gas 
collecting  tube  should  not  be  too  low.  Otherwise,  it  must  be 
pressurized  by  a  pump  in  order  to  supply  gas  to  the  pre-cooling 
system.  If  the  maximum  pressure  of  the  collecting  vessel  is  2 

atm  and  the  minimum  pressure  is  0.3  atm,  then  the  total  volume  of 

3 

the  vessel  must  be  2.7  -  3.6  m  .  If  low  temperature  air  is 
directly  exhausted  into  the  atmosphere,  then  the  temperature 
should  be  raised  to  protect  the  environment. 
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Figure  5.  A  2.4m  x  2.4m  High  Re  Number  Cryogenic  Transonic 
Wind  Tunnel 


1 .  atmosphere 

2.  five  Model  22  air  compressors 

3.  air  pre-cooling  systems  (two  in  parallel) 

4.  high  pressure  gas  storage  vessel 

5.  temperature  stabilizing  device  (636  tons  of  LF5 
corrosion  resistant  aluminum  alloy) 

6.  acoustic  damper  valve 

7.  pressure  stabilizing  box  150  atm 

8.  hypersonic  wind  tunnel 

9.  safety  and  protection  device 

10.  atmosphere 

1 1 .  shock  wave  2 

12.  atm  exhaust  gas  collecting  vessel  (30,000ni  ) 

13.  wave  stabilizer 

14.  isobaric  flow  throttle  in  the  stabilized  section 
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Abstract 

In  order  to  achieve  the  full  scale  Reynolds  number  of  aerocraft  in 
model  testing,  various  high-ReynoIds  transonic  wind  tunnels  are  being 
developed  abroad.  In  this  paper,  a  cryogenic  high-Reynolds  number  tran¬ 
sonic  wind  tunnel  with  pre-cooled  and  restricted  flows  is  presented. 
The  principle  that  air  temperature  falls  down  through  a  flow  restrictor, 
which  is  also  a  regulator,  is  applied.  Air  from  compressor  is  first  cooled 
to  215 "K,  and  then  enters  into  the  pressure  vessels.  During  the  wind  tun¬ 
nel  operation,  the  regulating  valve  must  be  controlled,  so  that  fluid  pres¬ 
sure  is  5atm  and  its  temperature  rs  15i°K.  Under  the  different  Mach 
number  condition,  the  different  temperature  and  pressure  may  be  used  to 
achieve  a  Reynolds  number  as  high  as  16.7x10’.  In  this  paper,  the  cool- 
system  of  the  wind  tunnel,  the  tunnel  operating  principles  are  described 
in  detail,  as  well  as  the  2.4 a  transonic  wind  tunnel  scheme. 


